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The  first  part  of  this  report  (I)  investigates  a  given  over¬ 
determined  set  of  lst-order,  partial  differential  equations  and 
reduces  them  to  the  proper  one  for  given  initial  and/or  boundary 
conditions,  in  a  natural  representation  for  finite-difference 
descretizations  and  calculations.  Alternative  variational  formu¬ 
lations  are  then  developed  and  the  same  questions  (of  determinancy 
and  well-posedness)  asked  with  respect  to  finite-elements 
descretizations  based  on  them. 

We  then  move  (Part. II)  to  the  actual  solution  of  the  Tricomi 
problem,  based  on  the  variational  formulations  proposed  and  on  finite 
elements  descretizations.  The  problem  involves  a  mixed-type 
(elliptic-hyperbolic)  domain,  and  serves  as  a  model  for  transonic 
flows,  when  represented  in  the  hodograph  plane. 

The  following  3  parts  (III,  IV,  V)  are  related  to  ai  new  design 
procedure,  proposed  for  transonic  airfoils  and  aircraft.  The 
solution  of  a  non-standard  elliptic  boundary-value  problem  is 
involved  and  its  mathematical  (III,  IV)  and  numerical  (V)  aspects 
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are  investigated.  The  latter  part,  exposing  and  implementing 
the  new  procedure,  is  preliminary.  It  needs  further  development, 
yet  may  hold  versatility  and  possibilities  beyond  other  design 
procedures  currently  available. 


REDUNDANT  SYSTEM  AND  VARIATIONAL  PRINCIPLES 


IN  CONTINUUM  MECHANICS 


Nima  Geffen 


Abstract 

NJ  The  mathematical  modeling  of  a  physical  continuum  by  a  set  of 
first-order  partial  differential  equations  is  commonly  redundant, 
i.e.,  it  includes  more  equations  than  unknowns.  A  simple  analysis 
gives  the  reduced  set,  which  fully  determines  the  field  quantities 
satisfying  appropriately  prescribed  initial  and  boundary  conditions. 
An  algebraic  look  at  the  system  provides  the  appropriate  choice 
for  descretization,  which  is  a  necessary  condition  for  solvability 
and  stability. 

The  same  is  tried  for  a  number  of  variational  representations 
of  the  field  in  terms  of  primitive  variables.  A  few  of  the  latter 
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1*  Forward. 

Two  forms  for  mathematical  codeling  of  physical  contic.ua  are 
daacrlbed  and  discussed  in  two  parts:  I.,  a  system  of  first  order , 
partial  differential  equations  for  the  field  components  and  quantities. 
II.,  an  integral  representation  of  the  fields  via  a  variational  prin¬ 
ciple  j  a  functional  of  the  field  quantities  becomes  stationary  at  the 
solution -point  of  the  differential  system  which  constitutes  it's 
Kuler-Lagrange  conditions.  The  two  approaches  are  common  in  physics 
(e.g.  the  Newtonian  and  Lagrangian  representations  of  mechanics) 
Hamilton's  principle  of  least  action)  and  both  have  been  extremely 
useful.  Illuminating ,and  necessary  in  many  and  varied  areas,  both  in 
classical  (e.g.  mechanics,  electrodynamics,  relativity  e.g.  [1],  [2], 
[3])  and  quantum  applications  [1], 

The  common  presentation  (e.g.  [1],  [2],  [3],  [4])  in  theoretical 
physics  ignors  the  question  of  well  posedness.  The  principle  of 
least  action  in  mechanics  is  formulated  as  a  boundary  value  problem, 
where  the  initial  and  final  locations  and  times  are  prescribed 
while  in  fact  it  is  an  initial  value  problem  where  location,  time 
and  velocity  are  given  initially  and  evolve  according  to  the  laws 
of  motion.  The  same  is  true  for  the  standard  variational  representa¬ 
tion  of  the  electromagnetic  field,,  relativity  and  quantum  electro¬ 
dynamic  . 

For  analysis  and  simulation  of  a  particular  physical  or  engin¬ 
eering  situation, appropriate  initial  and  boundary  conditions  have  to 
be  considered,  the  problem  has  to  be  well  posed,  i.e.  to  be  solvable 
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uniquely  end  steely, end  the  model  has  to  admit  -tome tries  and  non* 
smoothness  of  interest  in  practice. 

Initial  and  boundary  conditions  arc  Incorporated  In  the  follow* 
ing  analysis  for  both  the  differential  and  variational  formulations. 
A  few  of  the  variational  formulations,  we  believe,  are  new,  as  well 
as  viewing  the  redundancy  in  their  Euler  equations  as  a  possible 
cause  for  computational  instabilities.  The  initial  motivation  for 
this  work  came  from  finite -difference  [5]  end  finite-element  [6] 

simulation  of  problems  in  transonic  aerodynamics. 

The  question  and  system  treated  are  elementary  and  so  general, 
the  analysis  and  answer  so  simple,  to  be  judged  trivial,  if  not 
for  the  fact  that  the  question  does  come  up  occasionally  and  the 
answer  not  always  immediate.  Essentially  the  same  problem,  we  be- 
leive, has  been  treated  recently  (and  has  come  to  our  attention  wn.'ie 
writing  this  note),  in  a  completely  different  context  for  electro¬ 
magnetic  field  equations  ([  7],  [8]).  It  also  seems  to  be  re¬ 
lated  to  other  recent  approaches  [93  and  to  variational  formulation 
and  finite  element  analyses  of  analogous  continuum  systems  of  en¬ 
gineering  importance  [10 ],[!!],[  12] . 
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C.  General  Equations. 

Continuum  mechanics,  electrodynamics  and  other  physical  theories 
can  be  modelled  by  various  specialisation  of  the  following  set  of  first 
order,  partial  differential  equations  specifying  the  source  distri¬ 
bution  of  the  vector  field  A: 

U)  V.A  -  0 

and  the  curl  of  the  vector  field 

(2)  VXU-W 

for  n  +  1  independent  variables: 

U)  X  -  (x0,x1>(..,xn) 

where  may  designate  time,  (or  a  time-like  coord¬ 
inate). 

a  dependent  variables: 

(ii)  u  *  u(x)  »  (u1,...,ufl)(x0,x1,...,xn) 

(ill)  A  =  A( x,u)  =  (Aq,...,A^) 

(  iv)  G  =  G(x,u) 

(v)  W  =  W(x)  =  (W1,...,Wn). 


and: 
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If  u  appears  linearly  la  (111),  i.e.: 
aA  &a 


(vi)  —  •  ■  cij(x)  *hare  the  c^(x)  Ao  not  depend  on 

u,  which,  in  addition, does  not  appear  In  (iv)  i.e.: 


0  »  Q(x)  the  systea  (1),  (2)  it  linear. 

If  0  •  G(  x*u)  the  systea  it  teal-linear.  If  u  appears 
In  A  nonlinear ly,  i.e.: 


(vil) 


dA  &A. 

»  •  ’  Vs**> 


the  syttea  it  quasllinear. 


The  vorticity  W  has  to  satisfy  a  compatibility  condition, 


(an) 


V.W  -  0 


and  not  all  vorticity  distributions  are  admissible. 

The  system  (l),  (2)  can  be  extended  to  coupled  fields 
u(x),  2<x): 


(2«) 


{V  X  U  -  w 
v  x  v  «  n 


(2a*)  V.W  -  0 
(2a")  V.Q  -  0 


where  «  A^(x,u,v),  e.g.  the  electromagnetic  field  equations. 


The  system  of  equations  (l),  (2)  is  quite  general,  it  can  be 
linear  and  nonlinear,  elliptic,  hyperbolic  or  mixed, with  smooth 
or  non  smooth  solutions.  The  independent  variables  x^  may  desig¬ 
nate  space  and  time  coordinates  and  different  kinds  of  initial  and/or 
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boundary  coalition*  any  be  appropriate  for  different  problem*.  Higher 
order  equation*  nay  be  put  into  this  fora,  application*  include  problem* 
from  electro-magnetic  field  theory,  fluidynaalca,  and  plaaaa-dynaaics, 
including  flows  with  shocks. 


Initial  and  Boundary  Conditions. 

A  unique  solution  to  a  field,  governed  by  equations  (1),  (2)  is 

depicted  for  appropriate  initial/boundary  conditions  ,  where  a 

numer  (possible  all)  of  the  field  components  u .  sire  given  on  parts 

N  1 

(possibly  all)  of  the  boundary  afi  =  S  2Q  of  the  region  Q 

3 

occupied  by  the  field,  i.e. 

(3)  Uj(x)  *  Ij(x)  on  1  <  £  <  N. 

The  specific  form  of  (3)  depends  on  the  problem  and  on  the  type  of 
the  partial  differential  system  (1),  (2)  (e.g.  pure  boundary  values 

for  time  independent  elliptic  problems,  initial  values  for  the 
Cauchy  hyperbolic  case). 


3.  Redundancy  and  Reduction 


The  system  (1),  (2),  includes  (n  +  l)  first  order  partial 
differential  equations  for  the  in  unknown  coup  one  nts  of  the 
dependent  vector  u.  The  n  equations  (2),  however,  are  not 
linearly  independent.  In  order  to  solve  a  problem,  the  system  has 
to  be  reduced  to  a  'posed'  system,  a  necessary  condition  that  is  the 
number  of  equations  equal  to  the  number  of  unknowns. 

(a)  Potential  formulation. 

For  smooth  solutions,  one  has  traditionally  defined 
auxialliary  functions:  the  scalar  and  vector  potentials,  which  for 
the  system  (1),  (2)  can  be  written  as: 

(4)  A  =  -VxHg(0 

(5)  u  =  VO  +  f(X) 

where  ^(u,x)  is  a  vector  potential,  $  is  a  scalar  potential  and: 
g  and  f  are  any  solutions  of: 

(6)  Vog=G,  Vxf  =  W 

Upon  substituting  (k)  and  (5)  in  (1),  (2)  one  gets  a  2nd  order  system  of  1 
equations  for  the  k  components  of  the  scalar  and  vector  potentials  (4>  ,A) 

The  redundancy  in  the  first  order  system  is  exchanged  for  a 
redundancy  of  a  different  kind  in  the  higher  order  system,  the 
latter  admitting  a  family  of  superfluous  solutions,  all  describing 
the  same  field.  This  is  embodied  in  the  so-called  gauge  invariance 
for  the  common  scalar  and  vector  potential  representation,  and  the 
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actual  fields  (unique  for  determinate  well  posed  problems)  are 
obtained  with  the  aid  of  additional  gauge  conditions  and  an  appro¬ 
priate  consideration  of  initial  and  boundary  conditions  if  needed. 

In  classical  applications,  where  the  system  is  deterministic 
and  all  field  quantities  can  be  measured  directly  the  field-potentiax 
correspondence  is  clear  and  unique.  The  potentials,  although  pro¬ 
viding  a  very  useful  and  elegant  framework  can  nontheless  be  consid¬ 
ered  an  auxiliary  mathematical  entities  and  the  gauge  transformations 
a  mathematical  artifact,  while  the  true  physical  information  is 
completely  contained  in  the  fields  themselves.  Potential  representation 
has  been  immensely  useful  in  theoretical  physics  and  continuum  mechanics 
however,  the  primitive  field  variables  are  sometimes  preferable  especial 
ly  for  rotational  fields  and  for  flows  with  shocks,  for  numerical  simu¬ 
lations  based  on  both  finite  differences  and  finite  elements  descretiza- 


tions . 
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(b)  Addition  of  a  new  variable. 

A  suggestion,  to.  equalize  the  number  of  equations  and  unknowns 
and  still  solve  directly  for  the  fields  is  to  add  a  new  dependent 
variable : 


v  =  V  •  u 

and  using  the  relation: 

V  x  V  x  u  =  •  u)  - 

replace  equation  (2)  by  its  rotor: 

(7)  -  V  x  W, 

which,  with  equation  (1)  gives  (n  +  2)  equations  for  the  (n  +  2) 
unknowns  (u,v) . 

The  formulation  above  has  been  suggested  (and  used)  by  M.  Mock 
for  computational  purposes,  with  a  staggered  mesh  for  (u,v)  (to 
avoid  decoupling  of  the  descretized  equations  for  u  and  v)  to 
solve  boundary  value  problems .  This  formulation  is  not  directly 


extendable  to  the  unsteady  case. 


Direct  field  formulation. 

The  auxiliary  function  formulations  invariably  rai~e  the  order 
of  the  equations  to  be  solved;  the  first-order  system  becomes  second 
order.  This  requires  a  higher  degree  of  smoothness  for  the  solution 
function  and  may  be  a  draw-back  for  numerical  analysis  and  calcula¬ 
tions.  Thus  although  many  large  computer  simulations  are  based  on 
'potential'  formulations, a  direct  solution  of  the  first  order  system 
has  been  found  beneficial  [5],  sometimes  essential  [9],  especially 
for  'initial'  rather  than  boundary  value  problems,  and  for  non¬ 
potential  flows  (e.g.  Euler  solvers  in  aerodynamics). 

The  question  is  how  to  choose  the  'right'  n  rotationality 
conditions  that  will  give  with  the  continuity  equation  (1),  the 
'right*  (n  +  1)  x  (n  +  1)  system.  A  simp'*'}  treatment  for  irrota- 
tional,  3-dimensional  fields  (worked  out  for  the  problem  in  [5]),  is 
®Lven  in  [13].  It  is  extended  here  for  the  general  case  (equation 
(1),  (2)),  where  the  field  may  have  sources  and  be  rotational  (e.g. 
an  electromagnetic  field  with  moving  charges,  flow  behind  a  curved 
shock,  motion  of  reacting  gases.) 
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Choice  of  Equations  for  ?  apace-Dimensions. 

Consider  the  case  of  3  dependent  and  3  Independent  variables 
.(e.g.  Cartesian  space  coordinates).^ 

We  get: 

X  -  (x^XgfXp  =•  (x,y,z) 

H  *  (VVV  *  (u»v»w)  (x,y,z) 

A  -  a[2\  a<3>)  (u,v,v;  x,y,2) 


The  system  (1),  (2)  contains  k  equations  for  the  three  unknown 
functions  (u,v,w),  a  redundant  system  for  a  well  defined  field. 

In  addition: 

(2a)  B  0t 

x  y  z 

*  This  treatment  can  be  directly  extended  to  time- dependent  and 

larger  systems.  For  example,  for  the  three  dimensional,  time  dependent 

case, equation  (1)  has  an  additional  term,  with  an  additional  equation 

* 

(e.g.  an  energy  equation) J  Initial  conditions  are  given  everywhere. 
Equations  (2)  remain  the  same  and  the  whole  treatment  carries  over. 


Homtr,  for  twice  continuously  differentiable  (u,v,w)  and  once 
continuously  differentiable  t/2\  W^)  we  gat  by  differentiat 
ir-  equating  mixed  derivatives,  and  substituting: 


^tui) .  vn .  Uyt .  «43) 


•ubstituting  u  from  ( ii) : 
z 


-  w*y  ■  *i2)  -  >43) 

(Tz  -  V*  -  ^  +  "i33 


using  ( 2a) : 


and  integrating  with  respect  to  x: 


wy  -  *z  s^1^(xly,z)  -  w^Cx^y.z) 


Thus  equation  (i)  is  implied  by  equations  (iii)  and  (ii) 
provided  that  is  given  at  x  *  Xq  for  all  (x,y). 

In  the  same  manner  each  two  of  the  equations  determine  the  third, 


for  which  appropriate  initial  conditions  have  to  be  given. 
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tr»>taent. 


Written  in  a  matrix  font,  aquations  (1),  (2)  are: 


(8) 


4X)' 

0  0  0 
0  0-1 


0  10 


’u' 

■ 

u 

V 

+ 

0  0  1 

V 

w 

V J 

X 

0  0  0 

-10  0 

r. 

f  (3)  (3)  (3)1 

Au  Av  Av 

0-10 

10  0 

0  0  0 


a 


3 

S  ax 

i-1  0*1 


(1) 


lA) 

^2) 

l|(3) 


or: 

(8a)  Uy  +  F 


with  the  corresponding  (4x3)  matrices  C  and  4  forcing  functions 

r. 

me  system  (8)  is  again  redundant,  and  one  of  the  last  3 
equations  has  to  be  deleted  to  reader  (4x3)  coefficient 

matrices  C  into  (3  x  3)  matrices  C.  The  choice  is  directed 
by  the  marching  ’’time  like"  direction  for  which  initial  conditions 
are  given,  for  which  the  ^  matrix  has  to  be  invertible,  hence 
nbnsingular,  hence  with  no  row  (or-columi$  of  ceros.  This  auto¬ 
matically  rules  out  one  equation:  when  (u,v,w)  are  given  at 


x  •  Xq  («.g.  Xq  -»  -•  for  steady  flow  about  an  obstacle),  tb«  first 
lrrotationality  condition  has  to  bs  knitted.  A  marching  procedure 
along  the  x  axis  requires  the  inversion  of  the  nota-redundant  matrix 
C(  x)  such  that 

(8b)  u,  -  *  ? \)  *  i t*h*z. 

In  the  same  manner,  integration  schemes  along  the  y  and  s 
directions  require  the  omission  of  (Si)  and  8iii),  respectively, 
as  a  necessary  condition  (not  always  sufficient!)  for  a  well  defined, 
stable  scheme,  (as  is  seen  in  (5l). 
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4.  Exaaplea . 


(a)  Rotational)  steady  aerodynamics . 

Substituting  for  p(u2)  in  (la)  and  eliminating  the  time 
derivatives  one  gets  (e.g.  [5]): 


(10)  (a2  - u2)ux  4-  ( a2  -  v2)vy  ♦  ( a2  -  v2)v;uw ( w  +  *2)  -  ^(uy  +  vx)  -  v*( 

(i)  y  -  U  -  (6) 

X  s 

(ii)  Vy-V^-l/2* 

(iii)  U  -  V  - 

y  * 


i.e.  4  equations  for  the  3  unknown  functions  (u,v,w). 

The  potential  formulation  has  rot  been  found  feasible  for 
rotati' nal  flo.'s  (except  for  the  two  dimensional  case,  where  the 
■vector  potential  is  the  stream  function),  and  the  first  order  system 
is  solved  directly  for  this  easel?]. 

The  matrix  fonnutaticn  for:  u  ■  (u,v,v)»  includes: 


-(x)  B 


J  a2  -  u2  *uv  -uw 
0  0  1 

0  0  0 

0-10 


and  similarly  CVJf;and  C 


(y). 


,(*) 


r  *  (O)V^1',  J>2\  w^). 


and  the  equation  containing  the  row  of  0's  in  the  marching  direction 
(for  which  "initial"  conditions  are  known)  is  zo  be  eliminated. 


w  +v  )= 
*♦  -*  7 
J  ** 
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( al)  Irrotational,  unsteady  aerodynamics 

We  set: 

*  -  (vwV  ■ 

u  «  (u^Ugjup  «  (u,v,v) 

A  -  ( p,  u) 

,  2s 

P  -  p(u  ) 

G  «  0,  V  *  0 

and  get: 

(la)  Pt  +  (pu)x  +  (pv)y  +  (pw)z  n  0 


(2a)  V  x  u  «  0 

o  *  p(u2). 


Deriving  u  from  a  potential  *  :  u  *  V  $  satisfied  (?a) 
identically.  In  addition,  *  has  to  satisfy  a  2nd  order  nonlinear 
equation  derived  from  (la). 

The  first  order  system,  in  the  matrix  form  (3a)  is: 


u  ■  (v,u,w)  ,  F  ■  (0,0, 0,0) 


,(x) 


0  U  *  0 

u 


0 

0 

0 


0 

0 

1 


'•'1 


0  i 
-1 
0 


(y) 


and  similarly  C  and 
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(t)  Maxwells  equations  in  vacuum. 

The  electromagnetic  field  equations  are  [1]: 

(11)  V  •  H  »  0  (13)  7x1.^ 

(12)  V  •  E  ■  4tT0  (14)  7x  H*  ;E,  ( £  *1*1) 

Where  E,  H;  p, designate  the  electric  and  magnetic  fields, 
the  electric  charge  and  electric  current  distributions,  respectively. 
They  are  of  the  for®  (la),  (2a)  where: 

^  -  H,  Ag  -  E , 

0X  -  0,  Qg  -  4*  p 

W,  -  -  E.  +  ~  j  ,  W0  -  -  -  H, 

X  C  “X  C  *  2  C  “t 

Equations  (U)  -  (14)  are  not  independent:  for  smooth  fields 
(twice  differentiable)  equations  (13)  implies  (V  •  H)^  =  0  and 
equation  (11)  holds  for  all  times  if  it  holds  for  t  =  tQ,  thus  it 
is  but  an  initial  condition  (at  most),  that  has  to  be  satisfied  by 
H  at  t  «=  0. 

The  electric  charge  density  and  currents  (p  and  respect¬ 
ively)  cannot  be  prescribed  arbitrarily,  since  equation  (14),  and 
(12)  imply  a  constraint  on  their  source  terms  p  and 

(15)  p^.  +  ^  V  •  ^  =  0  (continuity  of  charge) 

which  has  to  be  satisfied  at  all  times. 

From  the  accounting  (and  counting)  point  of  view  for  3  space 
dimensions,  the  system  includes  8  equations  for  unknowns 
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( three  components  for  each  E  tad  H).  Equations  (13)  and  (14)  can 
be  viewed  as  evolution  equations,  while  (11)  and  (12)  are  mere 
constraints,  augmented  by  the  compatibility  condition  (15)*  This 
classification  is  certainly  not  unique  and  not  always  obvious  (e.g. 
in  the  stationary  case,  where  the  time  derivations  are  absent  and  the 
two  pairs  look  equivalent) . 
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Reduction  to  a  minimal  set* 

( i)  Potential  formulation 

for  Maxwells  aquation a  tha  standard  analysis  is  done  via  tha 
scalar  and  vector  potentials:  ♦  and  A  respectively,  ( a.|.  [11]), 
where: 

1  ^ 

E  ■  -  r  rr  -  ™  »  h  «  v  x  a. 

Written  in  terns  of  ($  ,A)  Maxwell's  system  reduce  to  1  scalar 
equations  for  the  U  components  of  * '  and  A. 

The  resulting  equations  are  higher  order,  and  admit  a  wider 
family  of  solutions,  all  equivalent  under  gauge  transformations: 

» •  •  ♦  - 1 
A’  .  A  +  Vf 

or  in  4  dimensional  notation: 

(A,  -♦)  '  -  (A,  -♦  )  +  [7f,  -  “  ft3. 

The  potential  formulation  for  the  electromagnetic  field  is  endowed 
with  a  beautiful  structure,  lucid  transformation  properties  and 
striking  accessible  information  content  (e.g.  decoupling  into  2nd 
order  wave  equations  for  each  of  the  components)  unfolding  the 
wealth  of  electromagnetic  waves  reference  [1], 
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(ii)  Adding  new  variables 
Setting; 

(16)  h  «  V  •  H  ,  e  =  V  °  E 

to  be  added  to  the  system  (11)  -  (15)  gives  the  additional  depend¬ 
ent  variables  required,  and,  again,  a  second  order  system  obtained 
by  substituting  the  expressions  (15)  in  (13)  and  (1^).  The 
system  (11)  -  ( 16)  the  sane  number  of  equations  and  unknowns,  at  the 
cost  of  2  additional  dependent  variables  and  a  2nd  order  equations. 

(lii)  Matrix  formulation 

Equations  (ll)  -  (14)  when  spelled  out  (in  Cartesian  coord¬ 
inates,  say)  can  be  readily  put  in  a  matrix  form  (J+),(5).  The 
criterion  of  a  non-singular  coefficient  matrix  in  the  marching 
direction  applies,  resulting  in  an  appropriate  elimination  of  the 
redundant  equations.  No  new  variables  are  added,  and  the  system 
remains  first  order,  hence  admitting  non-smooth  solutions  and  the 
treatment  of  nonlinear1  problems.  Combined  electromagnetic-fluid - 
dynamic  flows,  as  well  as  other  field  combinations  can  be  formulated 
reduced  and  discretized  accordingly. 


6J(  V)  a  0  <=  V  =  U  . 

For  well  posed  problems,  for  which  u(x)  is  unique 

5J(  v)  =  0  <==>  v  =  u. 

Variational  formulations  are  scalar,  short,  additive  (the  func¬ 
tionals  for  complex  systems  are  direct  sums  of  their  simpler  parts), 
invariant  under  appropriate  classes  of  transformations  and  are  often 
convenient  for  theoretical  analysis  and  for  numerical  simulations,  e.g. 
by  finite  differences  and  finite  elements  [lo],  *11 3,  [12];  the  integrals 
can  easily  be  descritized  and  approximated,  and  smoothness  require¬ 


ments  on  the  functions  are  less  stringent  than  for  the  corresponding 
differential  system.  This  last  point  is  most  important  from  the  numeri¬ 
cal  view  point  and  for  the  treatment  of  shocits. 

The  case  G  =  0  is  described  in  [l&’J.  The  functional  J(v,x): 


Vux  A  *  0 . 


A»VL 

—  u 


resulting  In: 


A  -  -VxX  . 


The  variation  is  done  on  all  v  for  which  J  is  defined,  coercive 
boundary  conditions  are  satisfied  on  the  boundary  =  oO  -  >  or 

for  which  ^||v  or  v||da  on  an  f 


For  the  non-sourceless  (or  sourceful)  case,  the  following  variational 
statements  hold: 

Theorem  1.  The  functional: 

(25)  J(v,>0  =  /  [L  -g  ♦  v  +  •  (Vxv-W)]  •  dx  +  /  ?ixv  •  dcr 

is  stationary  for  v  =  u  satisfying  (1),  (2),  (3)  provided  that 

(26)  VuxA-0 


where : 


V  •  g  =  G  or  g  =  V  G 


A  *  Vx^  +  g 


v  is  allowed  to  vary  over  all  functions  for  which  (25)  is  defined 
and  finite,  and  which  satisfy  the  coercive  B.C.  (3).  In  this  formu¬ 
lation  v  is  required  to  be  differentiable  while  _>  can  be  Just 


integrable , 


;-l  i  ,  .Hail-iili,  f  ■_  tiitikki-Lii  ,u 
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An  alternative  formulation  ia  obtained  by  integrating  the  second 
term  by  parts,  using  .the  vector  identity: 

v  *  (*xv)  «  V  *  7x>  -  >  •  Vxv 

substituting  in  (25)  ^  •  77xv,  we  obtain: 


Theorem  2.  The  functional 


(30) 


J(2>£) 


v  +  v 


^  *  W)dx 


is  stationary  for  v  «=  u  satisfying  (1),  (2). 

In  the  variational  formulation  (30),  the  Lagrange  multiplier 
T\  is  required  to  have  integrable  first  derivatives,  which  appear 
explicitly  in  J,  while  v  can  just  be  integrable. 

The'surfaci  term  drops  out,  and  the  solution  v  =  u  satisfies 
the  natural  boundary  conditions: 


A  xu  •  dcr  *  0  on  ^Q,  i.e.  » 


either 


*  ||  u  or  u  [|  i.e.,  u!!  fi  ^fl. 
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Examples: 

a)  Steady  electromagnetic  field. 

Maxwells'  equations  for  a  steady  state  can  be  written  as  [3]: 


1) 

V  •  E 

*  P 

ii) 

VxE 

.  0 

iii) 

V  •  B 

«  0 

iv) 

VxB 

-  A 

A  variational  statement  for  i) ,  ii)  is: 

(25')  J(E,*E)  *  f[(E2/2-£  •  E)  +  ^  xE  *  do; 

for  E  arbitrary,  or: 


(30’)  J(E)  *  /  (E^ /2  -  g-E)dx  for  E  irrotational 

JQ 

where:  g  is  any  solution  to 


V  •  g  ■  o  or  g  *  • 


The  corresponding  variational  foraulations  for  the  magnetic 
equations  iii),  iv)  are: 

J(B,>B)  -  /C32/2  +  >{B\vxB-<i)]dx  +  B^xB  •  da 


(25”) 
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b.  Steady  fluidynamics. 
The  differential  equations  are: 
V-(oa)  =  0 

0  *  pU2) 

Vx£  *  w 


and  the  corresponding  Lagrangian  L  and  functionals  are: 


(#) 


„  ik 

au 


pu  -  “r  or:  pu 


_  iL 

i  3U, 


J(u)  =  /  [  L  -*  X  •  (7xu  -  w)]  +  /  X  x  u  •  dcr 

.  Jo  '  JxT  ~ 


or 


(51) 


r«. 


J(u)  =  /  L-X«w  +  u*VxX 

J  n  “  " 


for  the  appropriate  function  spaces,  with  the  required  smoothness 
properties  and  constraints,  in  the  region  and/or  on  the  boundary. 

The  results  (25)  have  beer,  described  and  analyzed  in  [16],  [17] 
where  applications  to  specific  flows  are  given .  While  ( 2^) 
requires  that  u  be  smoother  than  X  ,  the  opposite  holds  for  (31), 
which  is  more  in  accord  with  the  physics  of  the  problem,  where  the 
field  variables  u  are  related  to  the  derivatives  of  the  'vector  pot¬ 


ential' 
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Redundancy  in  Variational  Formula t ions . 

Variational  formulation  have  been  widely  used  for  approximating 
field  solutions,  notably  by  the  Rita  and  finite  elements  methods.  The 
tacit  assumption  is  that  the  variational  problem  is  well  posed,  for  which 
a  necessary  condition  is  that  so  is  its  system  of  Euler  equations  and 
initial  boundary  conditions}  a  necessary  condition  in  turn,  for  stable 
solvability  is  that  it  be  non-redundant,  as  shown  in  part  I  above. 

The  variational  formulations  Eq.  (25),  (30)  are  reduadai  t  in  the 
sense  t>at  their  Euler  Equations  (1),  (2)  are  redundant.  Since  one  of 
the  rotational ty  conditions  (2)  is  implied  by  the  rest,  its  correspond¬ 
ing  j\  component  is  completely  redundant.  This  can  also  be  inferred  from 
Eq.  (2S)  which  is  analogous  to  (2)  (with  u  replaced  by  A  and  w  by 
(A  -  g)),  where  the  system  (29)  is  not  independent,  and  one  of  the 
equations  results  from  the  other  ones. 

This  redundancy  in  the  variational  formulations  is  less  apparent 
than  in  the  differential,  system  (where  a  simple  count  of  equations  and 
unknowns  signals  a  red  light  -  to  stop  and  look).  In  the  usual  approxi¬ 
mations  (e.g.,  by  the  Ritz  or  finite  element  method)  the  number  of 
equations  is  made  to  equal  the  number  of  unknowns,  obscuring  the  fact 
that  the  system  may  have  no  solution  ( for  incompatible  w  not  satis¬ 
fying  Eq.  (2a)).  Simple  round-off  Inaccuracies  may  turn  various 
procedures  unstable.  This,  we  believe  may  be  one  of  the  causes  for 
instability  experienced  in  various  mixed  and  hybrid  finite  elements 
procedures  (e.g.  [10]  -  [14]), 
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Two  rowdies  are  suggested: 

(i)  to  remove  one  rotationality  condition)  the  one  in  the 
time -like  direction  (or  mostly  so  for  nonlinear  problems) 

(ii)  to  add  another  term  to  the  functional: 

8V  •  v 

to  enforce  the  compatibility  of  the  w  components  and 
correct  any  deviation  thereof,  while  (i)  removes  one  of 
the  unknowns  (but  demands  attention  in  choosing  the  correct 
one  to  eliminate)  (ii)  adds  an  unknown  -  another  Lagrange 
multiplier  which  is  a  draw  back.  The  scheme,  however,  may 
be  used  as  a  check  on  simple  bench-mark  problems. 
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Concluding  Remarks. 

The  variational  formulations  in  primitive  variables  hold  for  a  vide 
class  of  multi-dimensional  nonlinear  systems  and  are  convenient  to 
descretize  for  arbitrarily  complex  geometries.  They  also  hold  promise 
for  treating  multi-dimensional  shocks  as  actual  jumps  with  correct 
jump  conditions  (which  most  other  methods  cannot  do),  by  treating  the 
integrand  in  the  functional  as  piece-vise  continuous  with  an  appropriate 
generalization  of  the  Erdmann  conditions  at  the  jump  front.  This  as  well 
as  further  stability  consideration  for  different  variational  problems  and 
the  relation  to  other  variational  formulations  in  primitive  variables 
(e.g.  [19l)  are  to  be  treated  in  a  following  report. 
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APH3NDIX 

The  simplest  redundant  system. 

The  redundancy  in  the  rotational! ty  equations  (2)  seems  somewhat 
analogous  to  the  trivially  simple  case  of  3  equations  with  2  un¬ 
knowns: 


L1)  a-jj*  +  aL2y  =  b1 
(L)  L2)  a21x  +  a22y  =  bg 


L3)  a^x  +  b^y  =  b^ 


The  first  2 

equations  with  2 

unknowns 

LI) 

»liX  +  a12y  *  *1 

12) 

&2]_x  a22^  *"*  ^2 

zaay  have: 

a)  a  unique  solution  (x,y)  if 

linearly  independent,  i.e., 

and  only 

(Aa) 

aU  a12 

a21  a22 

*  o. 

b)  a  one  parameter  family  of  solutions  (x(s),y(s))  if  they  are 
linearly  dependent,  i.e. 


(  a21*a22*^l'  =  A  ^2JL> 


which  means  that  one  equation  is  a  simple  multiple  of  the  other,  and 

y  =  .!iix  +  ^,.!2i!t  +  !2_ 

a12  al2  a22  a22 


which  may  be  parametrized: 


(*»y)  =  (  s j as  +  b)j  a  *  -  - —  =  - 


a12  a22 


c)  no  solution  if 


all  al2 
a21  a22 


b.  a. 


b2  S22 


Geometrically  the  equations  models  lines  in  the  plane,  and  the  3  cases 


u  - -  -  - - 

»  - 
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The  case  (a)  is  the  usual  (or  more  prevalent  one)  while  (c)  and 
(b)  are  exceptional. 

With  the  added  equation  the  (3x2)  system  (L)  the  system 
may  have 

A)  a  unique  solution  if 


and 


b. 

a,  _  b 

a _ 

1 

12  2 

22 

b2 

a22  b3 

a32 

a12  a21 

a22 

a21 

a22  a31 

a32 

*11  bl  a21  b2  j 

*21  b2  *31  b3  I 

all  a12  a21  a22 

a21  S22  a31  a32 


B)  a  one  parameter  family  of  solutions  if  all  three  equations  are 
linearly  interdependent,  or  any  2x2  sub -determinants  for  the  4  x  3 
matrix  (a._,,b.)  is  zero,  i.e.,  at  least  2  of  the  equations  are 
multiples  of  each  other. 

C)  no  solution,  when  none  of  the  conditions  above  holds. 


. .uirLJitLiL*bJui^i . .....  1  J- *  j jA..,  ,  ,  .  . ..  ■:,;J,-.;-.ili  i.i  .rii 
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A:  the  three  lines  L^,Lp>L^  intersect  at  one  point. 

A':  two  lines  (1^,1^)  coalesce,  intersecting  at  a  point. 

B:  three  coalscing  lines  L  :  (x,y)( s) . 

C:  three  lines  intersecting  at  3  points  -  the  most  prevalent 
situation  of  no  solution  to  the  system. 

The  case  discussed  in  the  paper  Is  analogous  to  the  case  A,  where 
presumably  there  exists  a  unique  solution,  satisfying  the  too  many,  but 
congruent  equations  "meeting  at  one  point".  The  problem  is  that  due 
to  inaccuracies,  the  equations  are  not  exact,  which  immediately  shifts 
us  to  the  situation  C  of  no  solution. 

This  simple  model  of  redundancy  readily  expends  to  larger  "rectangu¬ 
lar'  (rather  than  'square')  systems,  with  their  geometrical  representa¬ 


tion. 
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The  4  linear  equations  with  3  unknowns: 

(P)  ‘yXj  »  i  =  l,...,4,  J  =  1,2,3. 

describe  4  planes  (P^,...,P^)  in  3  dimensional  space. 

Depending  on  the  (4x3)  matrix  of  coefficients  a^j  and  the  vector 
b^,  or,  equivalently,  the  matrix  the  Planes  may: 

I.  Coalesce  to  one  plane  (if  all  the  equations  are  linearly 
dependent) . 

II.  Intersect  along  1  to  4  lines 

III.  Meet  at  one  point,  vhich  forms  a  unique  solution  [x^]  , 
i  »  1,2,3- 

This  last  case  is  the  most  special  and  requires  special  compatibility  con¬ 
ditions,  which  are  violated  by  round-off  errors  even  for  supposedly 
compatible  systems. 


SADDLE  POINT,  APPROXIMATION  FOR  THE 


TRICOMI  EQUATION 

by 

Sara  Yaniv,  Frieda  Loinger  and  Nima  Geffen 


Abstract 

The  Tricomi  equation  is  solved  in  a  mixed  elliptic-hyperbolic  domain 
by  a  finite-element  implementation  of  a  new  variational  principle.  The 
nonuniformly  elliptic  part  is  solved  by  using  a  projection  of  the  boundary 
condition  in  the  hyperbolic  region  onto  the  parabolic  line.  This  leads  to 
a  convergent  solution  which  is  extended  into  the  hyperbolic  part  with 
characteristic  finite  elements  by  solving  either  a  Goursat  or  a  Cauchy 
boundary  value  problem. 
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1.  INTRODUCTION 

The  Tricomi  equation,  one  of  the  simplest  equations  of  mixed  ellip¬ 
tic-hyperbolic  type  and  of  interest  in  physical  applications,  representing 
small  perturbations  transonic  flows  in  the  hndograph  plane,  has  been  inves¬ 
tigated  extensively,  e.g.  [1],  The  Tricomi  problem,  where  boundary  data 
are  given  on  an  elliptic  arc  and  one  characteristic  in  the  hyperbolic  half 
plane,  can  be  separated  into  a  nonuniformly  elliptic  problem,  where  bound¬ 
ary  data  are  given  on  the  elliptic  arc  and  the  parabolic  segment  'closing' 
it,  and  a  nonuniformly  hyperbolic  triangle  bounded  by  a  parabolic  "bai is" 
and  two  characteristics  emanating  from  its  sides  (figure  1).  Other  bound¬ 
ary  conditions  have  also  been  treated  [1], 

An  analysis  of  a  complete  linear  boundary  value  problem  of  whatever 

type  and  boundary  conditions  is  given  in  [2],  with  conditions  for  well 

losedness  and  a  weak  formulation.  The  last  analysis  is  constructive  and 

directly  amenable  to  discretization  of  a  mixed  elliptic-hyperbolic  domain, 

using  the  same  procedure  in  both  parts,  in  a  Galerkin-type  weak  formulation. 

A  finite-difference  scheme  based  on  this  weak  formulation  is  used  in  [6] 

]/-> 

with  convergence  0(h  “)  and  triangular  linear  finite  elements,  with 

2 

conve  -ce  0(h“)  in  [7],  where  a  comprehensive  survey  of  methods  and 
an  extensive  bibliography  are  also  given.  The  convergence  reported  in 
[7]  is  better  than  the  theoretical  estimates  in  [3]  for  the  same  problem. 


-  43  - 


A  direct  finite-element  analysis  and  estimates  of  the  nonuniformly  ellip¬ 
tic  part  is  given  in  [4].  Finite-elements  calculations  using  patched 
variational  principles  are  given  in  [12]  with  "experimental"  second  order 
accuracy  for  the  tested  problems.  Finite  difference  calculations  have 
also  been  tried,  e.g.  [S], 

Another  mixed  type  equation  solved  recently  is  the  Lavrent'ev- 
Bitsadte  equation  [14],  where  the  boundary  condition  on  the  hyperbolic 
region  are  transferred  to  the  parabolic  segment  and  the  equation  in  the 
elliptic  region  is  solved  by  the  finite  element  method. 

In  (IS)  a  weak  formulation  for  the  problem  based  on  different 
spaces  of  test  and  trial  functions  is  used  to  prove  existence, 
uniqueness  and  uniform  stability  of  the  approximate  solution. 
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2.  A  VARIATIONAL  FORMULATION 
Tricomi's  equatioit: 

*xx  -  *yy  ■  0  <» 

is  elliptic  in  the  lower  half  plane  y  <  0  and  hyperbolic  in  the  upper 
one,  y  >  o  (figure  1),  where  it's  real  characteristics  are: 


e  ■ 


♦  f  >'V2 


(2a) 


(2b) 


Writing  equation  (1)  as  a  first  order  system,  one  gets: 


(u.v)  =  (4>y,<j>y) 

(yu)x  -  vv  *  0  (la) 

uy  '  vx  =  0  *  (lb) 

We  now  look  for  a  solution  in  a  mixed  domain  fi  *  f.j  U  0,  where 
is  in  the  lower  (y<0)  elliptic  half  plane  with  its  boundary 
and  in  the  upper  (y>0)  hyperbolic  one  with  it's  boundary  Fj  U  F, 

where  1^  and  r,  are  two  characteristics  with  boundary  conditions  given 
along  and  one  of  the  characteristics,  say  Fj  (figure  1).* 

Other  domain  shapes  with  different  conditions  can  be  similarly  considered. 
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2.1  A  variational  formulation  in  the  elliptic  region. 
First,  let  us  consider  the  Dirichlet  problem: 


■  *yy  ■  f(x'y) 


in  a  domain  n  (v<0)  with  given  boundary  conditions: 


♦(x,y) 


Assuming 


U  ■  4 


v  *  4. 


then  the  operator  (uy)x  -  v  is  a  potential  [13,  page  35], 


The  variational  formulation  for  the  problem  is:  find  (u,v;X)  so  that 

r 

J(u,v; X)  =  |  [L(x,y,u,v)  ♦  X(u  -v  )]dxdy  (41 

>  (?'  •  * 

r  • 

-  21  u»F(x,v)dxdy 

JoJ 


where 


fX 

F(x,y)  ■  j  f(C,yldC, 


is  stationary,  for  all  functions  (u,v)  €  U*V  and  X  6  K. 


— tideu ‘.m iu-JJ.  ui,.  siflitotL. 
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^(x,y,u, v) 


is  the  Lagrangian  of  the  problem,  i.e. 


L  =  vu 
u 


L  =  -v 

V 


«  is  a  Lagrange  multiplier  which  is  the  stream  function  of  the 

problem  [13,  page  39]. 

The  variation  of  J(u,v;A)  for  fixed  A  giv« 
formulation : 


gives  the  following  weak 


r  r 


6J(u,v;A)  »  J  Jyudu  -  v«v  ♦  \(«uy-«vx)]dxdy 


C5) 


Adding  the  equation 


[  [iu*F(x,y)dxdy  *  o 


in  weak  formulation: 


'  vx  =  0 


Lta- 


vx)dxdy  »  0  for  VqCW, 


we  get  a  saddle-point  problem: 


t  * 

u 


[yuiij  w1  +  A  i  -  3x')  ]dxdy  =  f  jUjF  dxdy  (6a) 


v(u] ,Vj ;  €  UxV 


Jqv  *  Tx-* dxdy  *  0  Vq€W 


In  the  lower  half  plane  yu2  -  v2  *  0, 


hence 


U*V  •  t (u.vD/f  f [-yu2+v2*(u  -v  )2)dxdy  <  •, 

JnJ  y  x 


udx  +  vdy  ■  0}  , 
30 


l(u,v)| 


U*V 


[-yu2+v2* (u  -v  ) 2 ] dxdy 

&  y  x 


and,  since  a  depends  on  an  arbitrary  constant : 


W  -  {q  j 


r  r  2 

|q  dxdy  <  •  ,  I  I  qdxdy  ■  0  > 

nJ  V 


2  f  f  2 

191“  B  |q  dxdy 

JnJ 


This  formulation  gives  a  saddle-point  problem  [9],  which  is 
determined  as  follows: 

Given  f  €  V’,  g  €  r,  find  u  €  V,  <|»  €  W  such  that: 
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a(u,v)  ♦  b(v,iii)  *  f(v)  W€V 


b(u,4)  *  gU) 


V*€W 


r' 


Bre22i's  theorem  states  the  following  conditions  for  existence 
and  uniqueness  of  the  solution: 

Suppose  a  and  b  are  bounded  and 


inf  sup  !a(u,v)j  $  y  >  0  (8a) 

u€z  v€z 

lulv»l  »vi-v»l 

z  =  {u€V/b(u,*)'»  0  V$€W) 


SUP  5  k  1*1  W  V*€W»  k  >  0  (8b1 

vev  “V“v  w 


then  the  solution  of  (5)  is  unique  and 


lUly  *  M»IW  «  C(  |f  |y  ,  ♦  ||g||W,). 


It  is  easy  to  verify  that  problem  (6a),  (6b)  satisfies  (8a), 
(8b),  hence  has  a  unique  solution  which  is  the  only  solution  of  (3a), 
(3b). 


l.-.r.-f'-'W”' 


\  .  ...  -  49  - .. 

In;,  order  to  approximate  the  solution  cf  problem  (6a)  *  (6b)  using 
the  saddle-point  weak  formulation  or  the  saddle-point  variational 
principle,  one  should  use  finite-dimensional  spaces  ^  c 
and  W'h  c  satisfying: 


inf  sup  |a(u,v) |  5  t  >  0  ,  t  independent  of  h  (9a1 

u£Zh  v£zh 
iiup  =1  II  v|iv=l 

sup  |a(u,v)|  >0  V  0  ^  u  6  Zh  (9bl 

V^h 

Zh  =  (u€Vh;  b(u,*)  =  0  V4>ewh } 


sup 

v€vh 


!  b  (v,  iji) 
l|V|lv 


5  i 


l#lW 


(9cl 


v^ew. 


£  >  0,  independent  of  h, 


The  approximated  problem  is: 


a(u^,  v)  +  b(v,tph)  =  f(v)  we\'h 


b  ( ,  d> )  =  g(a) 


(10. 


V<t>€K, 


segment . 
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22  A  variational  formulation  in  the  hyperbolic  region. 

The  smooth  solution  in  the  elliptic  region  provides  mil  the 
data  on  the  parabolic  segment  which  is  needed  to  solve  either  a  Goursat 
problem  or  a  Cauchy  problem  in  the  characteristic  triangle  bounded  by 
Tl  U  r2  and  the  parabolic  basis. 

(x)  Consider  the  Goursat  problem: 

yux  -  v  ■  0  in  n2  U2a) 

u  -  v  ■  0 

y  * 

with  the  initial  conditions: 

■  * 

u(x,o)  ■  UQ(x)  -1  i  X  i  1 

and  either  u  /y~" ♦  v  | r  ■  f(y)  or  an  equivalent  condition  for  A:  (12b) 

‘l 

a  |  r  *  *\rr p(y)  dy 

The  variational  formulation  for  the  problem  is: 

Find  (u,v;X)  €  U  *  V  *  W  satisfying  (12b)  for  X  having  fixed 
(.unknown)  values  along  Tj,  for  which  the  functional: 


J(u,v;X) 


2  2 

[yu  -v  +l(uy-vx)]dxdy  ♦ 


X (udx+vdy) 


(12c) 


riUr2 


/ 
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is  stationary. 

This  variational  formulation  in  the  hyperbolic  region  is  a  con¬ 
tinuation  of  the  saddle-point  formulation  used  in  the  elliptic  region. 
But,  in  the  hyperbolic  domain  (y  >  0)  conditions  (8a),  (8b)  are  not 
satisfied,-  and  we  loose  the  proof  for  the  existence  and  the  uniqueness  of 
the  solution  of  the  variational  problem. 

(ii)  Cauchy  problem  for  the  Tricomi  equation.. 

The  given  boundary  conditions  for  the  problem  are: 


u(x, 0)  =  uQ(x) 
v  (x,0)  =  vQ(x) 


•1  {  x  f  l1 


(1*0 


An  equivalent  variational  formulation  [13  chapter  31  is: 

Find  (u.vjA)  such  that  (u,v)  €  U  *  V  and  satisfying  (13a) 


and  A  €  W  with  fixed  (unknown)  values  along  U  for  which 
the  functional: 


is  stationary. 


53 


3.  EXTREMAL  PROPERTIES  OF  THE  VARIATIONAL  FUNCTIONAL 

Conditions  (8a),  (8b)  given  by  Brezzi  [9]  are  satisfied  only 
in  the  elliptic  region.  These  conditions  are  connected  with  the 
following: 

Given 

t 

J(u,v;A)  =  j 

for  any  given  A,  the  extremum  of  the  functional,  noted  by 
(Uq(A),  Vq(A))  satisfies  the  equations: 


[vu“-v“+ A  (uy-vx)  )d'-;dy 


2yu  -  Ay  *  0 

2v  -  A  =0 
x 

hence  vu  -  =  0. 

•  x  y 

For  hyperbolic  and  mixed  type  equations  (up(A),  V0(A))  is  a 

saddle  point  which  may  be  specified  in  the  hyperbolic  case  in  what 
follows:  The  projection  of  the  functional  on  the  U  space  and  on 
the  V  space  separately,  yields  an  extremum  in  each  of  these  cases; 
the  projection  on  the  U  space  gives: 

■>  f  f 

6"J (u° (v, A) )  =1  '  2y  dxdy  ,  for  y  >  0, 

U  JP.J 
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so  u8(v,A)  is  a  minimum  in  U,  and  for  the  space  V,  the  second 

variation  of  J(u,v;A)|  is 

IV 

o  r  [ 

S~J (v* (u, A)) !  *  -I  2dxdy 

jv  JaJ 

so  v°(u,A)  is  a  maximum  in  V. 

Hence,  for  the  hyperbolic  region  the  '•point*' 

(u° (A) , v° (A) )  is  a  min  max  J(u,v;A). 

u€U  v€V 


For  mixed  type  equations  the  situation  is  more  complicated. 
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4.  APPROXIMATION 

We  use  the  variational  principle  (4)  to  solve,  approximately,  the 
Tricomi  equation  in  a  mixed  domain. 

The  boundary  conditions  for  the  4  examples  solved  are  chosen  so  that 
the  exact  solution  for  the  potential  <|>(x»y)  is: 


II.  <J>(x,y)  *  3x2  +  y3 

3  t 

III.  <Kx,y)  *  x  +  xy“ 

IV.  <Hx,y)  *  2x3y  +  xy4. 


The  numerical  results  are  illustrated  in  Tables  (I-V)  for  the  first 
example  (the  values  of  the  rate  of  convergence  of  the  other  examples 
behave  almost  the  same) . 

i)  The  Dirichlet  problem  for  the  Tricomi  equation  in  the  ellip¬ 
tic  domain.  Given  the  problem: 

y4>  -  $  =0 

/rxx  Tyy 

♦  -  f(x,y) 

an 

The  variational  functional  is: 

J(u,v;A)  -  f  [yu2-v2+  (u  -v  Jdxdy 
JQJ  r 


for  (u,v)  €  U  x  V,  X  €  W 
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ft  ->  -i  *> 

U  x  v  =  { ( u ,  v )  / I  [-yu‘*v‘+(u  -v  )fc]dxdv  <  «, 
J  >  x 


udx  ♦  vdy1  =  df} 

I  an 


W  =  {q/ 


JoJ 


q“dxdy  <  » 


•  L 


qdxdy 


oj 


for  this  problem  conditions  (8a),  (8b)  are  fulfilled,  hence  the  prob¬ 
lem  has  a  unique  solution. 

Let  n  be  the  rectangle: 


n  =  ( (x,y)/  -1  $  x  $  1  ,  *1  .<  y  {  0} 


we  divide  n  into  rectangles  of  size  h  and  take  different  finite 

dimensional  space  for  the  trial  functions. 

a)  The  first  attempt  is  to  use  bilinear  trial  functions  for 

for  (u,v)  and  X,  for  these  space  Brezzi's  conditions  (9a),  (9b)  and 

(9c)  are  not  satisfied;  this  is  concluded  from  the  instability  of  the 

numerical  solution.  For  these  trial  functions,  central  approximation 

for  A  and  A  are  received  and  therefore  the  value  of  A  at  one  cor- 

x  y 

her  determines  its'  values  only  at  every  second  point  in  each  direc¬ 
tion.  To  couple  the  equations  for 


'1 

I  ^ 

II 

ii 

i  3 
!  ^ 


u 

■a 

1 

i 
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all  grid  points,  the  values  of  '*  at  the  4  points  of  the  comer 

element  have  to  be  predetermined  (e.g.  by  Taylor's  expansion  about 

the  (-1,-1)  point  and  the  values  of  v,v  along  the  boundary). 

The  L,,  error  of  the  functions  u,v,x  for  the  4  examples 

given  above  is  approximately  0(h“)  for  the  Dirichlet  problem  (tee  Table  la). 

The  experimental  rate  of  convergence  is  reduced  to  1,6  if  the 

projection  (11)  is  used  along  the  parabolic  segment . (see  Table  Ilia). 

b)  Other  finite  dimensional  spaces  for  U^,  and 

for  which  (9a) ,  (9b)  and  (9c)  are  satisfied,  and  tried  (Appendix  A). 

Let  us  use  the  same  finite  elements  discretisation  of  ft  as  in  a  I , 

approximating  u  and  v  by  bilinear  trial  functions  but  piece-wise 

constant  functions  for  \  (intuitively,  this  leaves  us  with  the  only 

r  r 


arbitrary  constant  of  the  problem  which  is 


qdxdy 


0). 


For  this  approximation  the  solution  converges  to  the  analytic 
solution  for  both  of  the  problems;  The  Dirichlet  problem  and  problem 
with  the  integral  condition  (11)  (without  adding  artificial  condi¬ 
tions  for  stability). 

The  rate  of  convergence  of  the  L,  error  for  u  and  v  is  about 
1.4  and  1  for  >  (see  Tables  lb  and  Illb), 


iij  Goursat  conditions  for  the  hyperbolic  region. 


The  smooth  solution  in  the  elliptic  part  provides  the  values 
of  u  on  the  parabolic  line  A  B  (y  «  0,  -1  #  x  $  1),  which,  together 
with  given  boundary  conditions  along  the  characteristic  Tj  pro¬ 
vides  the  data  needed  for  the  Goursat  problem  in  the  triangle  ABC 
(well  posed  for  the  continuous  case).  The  region  is  now  divided  into 
isoparametric  triangles  formed  by  the  characteristic  mesh  in 
(figure  2)  .  Linear  trial  functions  in  £  and  n  (the  characteris¬ 
tic  variables)  are  used: 


where 


a*n  ♦  b*£  ♦  c 
m  m 


e 

m 


v®  •  (ue,ve,Ae) 


in  the  triangular  isoparametric  elements. 

The  scheme  is  inconsistent  near  the  parabolic  line.  To  overcome 
this  difficulty,  linear  elements  in  x,y  coordinates  and  linear  trial 
functions  are  taken  for  the  first  row  of  elements. 

The  functional  (12c)  is  discretized  and  the  nodal  values 
(ue,ve,Ae)  sought  from  the  discrete  algebraic  conditions  for 
stationarity. 

The  numerical  solution  of  the  same  examples  used  for  the  elliptic 
region  for  accurate  boundary  conditions  is  convergent  and  the  experi¬ 
mental  rate  of  convergence  of  the  L_,  error  is  about  1.5.  (see  Table  II). 
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(iii)  Cauchy  conditions  for  the  hyperbolic  region  02. 

An  alternative  way  to  treat  the  hyperbolic  region  is  to  consider 
(u,v,X) (x,0)  on  AB  from  the  elliptic  solution  as  Cauchy  (or  initial) 
data  for  the  hyperbolic  triangle  ABC,  and  use  the  functional  (13b) 
for  each  isoparametric  triangular  element  separately.  This  proce¬ 
dure  also  yields  a  well  posed  problem  with  explicit  algebraic  linear 
equations. 

The  solution  of  this  approximation  has  a  smaller  L2  error  than 
the  solution  for  the  Goursat  problem,  the  rate  of  convergence  of 
the  Lj  error  is  about  1.7  (see  Table  II). 

Note;  All  implicit  schemes  for  the  Cauchy  problem  sews  to  be 
unstable  and  explicit  schemes  with  Goursat  boundary  conditions  do 
not  converge. 

(iv)  Approximation  for  the  mixed  domain. 

As  has  been  shown  in  the  previous  chapter,  a  solution  for  the 

mixed  domain  is  available. 

1 

The  first  step  is  to  solve  the  problem  in  the  elliptic  region 
using  the  integral* projection  (11).  The  values  of  u  and  v  and  X 
along  the  parabolic  segment  are  achieved  by  this  calculation  and 
may  be  used  as  data  in  order  to  get  a  solution  in  the  hyperbolic 
region  by  solving  either  a  Goursat  problem  or  a  Cauchy  problem. 

The  results  are  illustrated  in  Tables  Ilia  and  Illb. 
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If  X  is  bilinear  in  the  elliptic  region,  we  get  large  errors 
in  the  values  of  X,  especially  in  the  last  example.  On  the  para¬ 
bolic  line  X  is  used  as  initial  conditions  for  the  Cauchy  problem 
but  not  for  the  Goursat  problem.  Therefore,  the  rate  of  conver¬ 
gence  is  better  using  Goursat  boundary  conditions. 

On  the  other  hand,  if  piecewise  constant  function  is  taken 
for  X  in  the  elliptic  region,  the  Cauchy  conditions  give  better 
results  than  the  Goursat  conditions .  The  schemes  for  the  Cauchy 
problem  are  explicit;  hence,  the  procedure  is  better,  and  finer 
meshes  for  the  hyperbolic  region  are  available. 

The  disadvantage  of  the  process  described  before  is  that  the 
integral  of  Bitsadze  (11)  needs  much  time  and  must  be  recomputed 
for  different  boundary  conditions. 

In  examples  II  and  III  difficulties  in  the  convergence  of  the 


problem  by  adding  a  constant  to  $(x,y)  so  that  the  new  function 


is  zero  at  x  ■  -1.  For  a  general  function  integration  by  parts 


must  be  done,  which  costs  more  computation  time. 

Another  disadvantage  of  the  Bitsadze  relation  (11)  is  that  it 
can  be  used  only  for  the  Tricomi  equation.  For  another  mixed  problem 
another  projection  formula  has  to  be  found;  therefore,  we  look  for 
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a  numerical  process  which  connects  between  the  hyperbolic  and 
the  elliptic  region. 

The  hyperbolic  region  divided  into  isoptranetric  tri¬ 

angular  elements  as  described  before  (Figure  2).  We  define  a 
local  variational  formulation  for  each  element  separately: 

Ji,j(u,v,X)  ■  J  ^  J  [yu2-v‘+X(uy-vx)]  dxdy  + 


X(udx  ♦  vdy) 


where  and  are  the  two  characteristics  bounding  the 


isoparametric  element,  and  use  linear  trial  functions  as  before. 

If  u<°>  -  the  solution  on  y  »  0  is  not  known,  the  explicit 
scheme  enables  us  to  find  a  connection  between  and  the  boun¬ 


dary  condition  along  1^.  For  every  j  we  can  find  a  matrix 
,  so  U™  ■  where  is  the  solution  at  the 

j-th  row,  and  is  a  matrix  depending  only  on  the  size  of  the 
mesh  and  not  on  the  boundary  conditions. 

Substituting  the  boundary  conditions  along  1^,  we  get  one 
equation  for  every  point  on  the  parabolic  line.  The  two  other 
equations  for  the  three  unknowns  (u,v,X)  will  come  from  the 
schemes  for  the  elliptic  region.  After  solving  the  elliptic 
problem  with  the  numerical  projection,  we  get  immediately  the 
solution  in  the  hyperbolic  region  using  the  explicit  scheme. 
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In  the  solution  of  the  elliptic  region  oscillations  of  the  values 
of  u  along  the  parabolic  segment  are  observed.  The  process 
converges  in  spite  of  the  oscillations.  In  order  to  demonstrate 
this  phenomenon,  the  error  of  u  is  calculated  in  the  usual 

L-  norm  and  in  the  weighted  norm:  a .  > .  (Table  IV). 

2  /^n/yu“dxdy 

The  oscillations  can  be  removed  by  averaging  the  values  of  u 
along  the  parabolic  segment,  and  using  the  smoothed  values  as 
Dirichlet  conditions  for  the  problem,  the  L2  error  of  u  is 
reduced  and  the  L2  error  of  vv  and  X  is  unchanged  (Table  V) . 

For  X  bilinear  in  y  <_  0,  the  rate  of  convergence  in  the 
elliptic  region  is  about  1.3  for  the  first  example  but  only  0.5 
in  the  last  example;  contrary  to  the  case  where  X  is  piecewise 
constant,  there  the  rate  of  convergence  is  the  same  for  all 
boundary  conditions  considered.  Therefore,  it  is  not  worthwhile 
fo  use  the  numerical  projection  together  with  X  bilinear  in 
the  elliptic  region  (Appendix  E) , 

Good  results  with  X  bilinear  in  the  elliptic  region  are 
achieved  oily  by  using  the  integral  relation  (11)  and  the  implicit 
schemes  of  the  Goursat  problem  in  the  hyperbolic  region;  hence, 
much  computation  time  is  needed  in  this  case.  Because  of  all  these 
disadvantages,  it  is  preferable  to  use  X  piecewise  constant  for 
all  the  problems. 
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5.  SUMMARY,  CONCLUSIONS  AND  REMARKS 

Tricomi  boundary  value  problem  in  a  mixed  domain  is  solved  on  a 
finite-element  mesh,  using  a  variational  formulation  admitting  arbit¬ 
rary  initial  and/or  boundary  conditions.  Convergent  solutions  are  ob¬ 
tained,  by  first  solving  the  elliptic  part  using  either  an  integral  condi 
tion  or  a  numerical  projection  of  boundary  condition  along  the  para¬ 
bolic  line,  then  extending  the  solution  into  the  hyperbolic  part  on 
characteristic  finite  elements,  solving  either  a  Goursat  boundary  val¬ 
ue  problem,  or  a  Cauchy  initial  one. 

The  advantage  of  this  approach  over  others  (e.g.  the  one  using 
Friedrich's  ingenious,  yet  non-trivial  procedure  from  the  practical 
point  of  view  even  for  the  simplest  cases)  is  that  it  is  conceptually 
simple,  does  not  require  further  transformations,  and  can  be  extended 
to  nonlinear  problems  and  to  multi-dimensional  and  unsteady  situations, 
at  least  in  principle. 


-  64  - 


REFERENCES 

[1]  Bitsadze,  A. V . , "Equations  of  the  Mixed  Type",  Pergamon  Press, 

1964,  (p.  86). 

[2]  Friedrichs,  K.O.,  "Symmetric  Positive  Linear  Differential  Equa¬ 
tions,"  Comm.  Pure  Appl.  Math.  11  (1958)  pp.  33-418. 

[3]  Lesaint,  P.,  "Finite  Element  Methods  for  Symmetric  Hyperbolic 
Equations,"  Numer.  Math.  21,  (1973)  pp.  235-244. 

[4]  Trangenstein,  J.A.,  "A  Finite  Element  Method  for  the  Tricomi 
Problem  in  the  Elliptic  Region,"  SIAM  J.  Numer.  Anal.  vol.  14, 

No.  6,  December,  1977. 

[5]  Ogawa,  H. ,  "On  Difference  Methods  for  the  Solution  of  a  Tri¬ 
comi  Problem,"  Trans.  Amer.  Math.  Soc.,  100,  (1961)  pp.  404-424. 

[6]  Katsanis,  T.,  "Numerical  Solution  of  Symmetric  Positive  Differ- 
erential  Equations,"  Math.  Comp.  22  (1968)  pp.  763-783. 

[7]  Azia,  A.K.  and  Leventhal,  S.H.,  "Numerical  Solution  of  Partial 
Differential  Equations  of  Elliptic-Hyperbolic  Type,"  Numerical 
Solutions  of  Partial  Differential  Equations  III,  Academic  Press, 
1976,  pp.  55-88. 

[8]  Geffen,  N.  and  Yaniv,  S.,  "A  Note  on  the  Variational  Formula¬ 
tions  for  Quasi-Linear  Initial-Value  Problems,"  Journal  of 
Applied  Mathematics  and  Physics  (ZAMP)  vol.  2",  1976,  pp.  833-838. 


li  -MVI 


Jte, 


•rrra'vv 


»  /safewfaaaaiaii 


-  65  - 

[9 ]  Brezzi,  F.f  "On  the  Existence,  Uniqueness  and  Approximations  of 
Saddle-Point  Problems  Arising  from  Lagrangian  Multipliers," 
R.A.I.R.O.,  1974,  p.  2,  pp.  129-151. 

[10]  Geffen,  N.,  "Variational  Formulation  for  Constrained  Quali-Lin- 
ear  Vector  Systems,"  Quarterly  of  Applied  Mathematics  October, 

1977,  pp.  375-382. 

[11]  Geffen,  N. ,  "Variational  Formulations  for  Nonlinear  Wave  Propaga¬ 
tion  and  Unsteady  Transonic  Flow,"  Journal  of  Applied  Mathematics 
and  Physics  (ZAMP)  vol.  28,  1977. 

[12]  Fix,  G.,  Gurtin,  M. ,  "On  Patched  Variational  principles  with  applica¬ 
tions  to  elliptic  and  elliptic-hyperbolic  problems,"  Num.  Math, 
vol.  28,  1977,  pp.  259-271. 

[13]  Yaniv,  S.,  "Variational  Formulation  For  Formally  Symmetric  Prob¬ 
lems  and  It's  Application  to  the  Tricomi  Equation,"  Ph.D.  Thesis, 
Tel-Aviv  University,  Israel,  1978. 

[14]  Deacon,  G.,  Csher,  S. ,  "A  Finite  element  method  for  a  boundary 
value  problem  of  mixed  type,"  SIAM  J.  Numer.  Anal.  vol.  16,  No.  5, 
October  1979. 

[15]  Aziz,  A. K. ,  Schneider,  M. ,  Werschulz,  A.,  "A  Finite  Element 
Method  for  the  Tricomi  Problem."  Num.  Math.  Vol.  35,  1980, 
pp.  13-20. 


-  66  - 


Appendix  A:  Stability  proof  for  the  elliptic  region. 


To  prove  the  stability  of  the  numerical  procedure  for  the  elliptic  region 
y  <  0, using  bilinear  approximation  for  u  and  v  and  piecewise  constants  for  X, 
the  following  condition  has  to  hold  [9]: 


sup  J 

(u,v)€UhxVh  n 


/ 

V  (uv  -  V  jdxdy 
))  y  * 


••  u  h 


for  every  !€  H. 


where 


element,  and 


U^xVjj  ■  { (u,v)/  u  and  v  are  bilinear  functions  in  each  rectangular 


rf  2  2  2 

[-yu  +  v  +  (uy  -  vx)  ]dxdy  <  »,  *dx  +  vdyj 


with  the  norm 


1 1 Cu,v) | |2UhXvh  =  J  J  t-yu2  +  u2  +  (uy  -  vx)2]dxdy, 


wh  =  {y 1^  is  piecewise  constant  in  each  rectangle 


j  4*2dxdy  <  «,  Jjydxdy  =  0} 


with  the  norm 


»  r  « 


Jjy^ixdy. 


P1 


'IKS?*-*™?- 
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Consider  the  patch  ot  four  elements  shown  in  figure  3  (the  same  holds 
for  any  number  of  elements) . 


Let  us  show  that  for  each  element  E. 


where 

and 


V  =  =  constant 


5*  i e  o 


it  Is  possible  to  find  a  pair  of  functions  (u,v)  €  U^xV^  for  which 


(15) 


f  f (Uy  -vx)  dxdy  «  Vi 


The  algebric  system  of  (15)  for  (u.v.)  i  »  0,1,... 8  is: 

9  J 


-(Vui) 

V  +V- 
0  o 

2k 

2h 

-<VUS> 

.  v.+v 
+  3  o 

2K 

2h 

Uo+U5 

2k 

V0+V7 

+  IF" 

VU1 

VV7 

2k 

'  IF 

^2 


and 


and  the  solution  is: 


4 

Vi 

1=1 


*  ¥„ 


0  . 
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f  VU1 
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*3  +  f4 


u +u- 
o  b 

11c 


vo+V3 


Vo+V7 

IK 


Let  us  choose  UQ  =  v0  a  °*  t  *  then 


jf<“y  -vx)2<lxdy  *  t  J1rijC'i  •  V2  S 


If  the  region  is  devided  into  more 
hence  a  is  independent  on  the  mesh. 


elememts.the  constant  a  is  not  changed. 


Calcuiating  \ i Cu»v) 1 1  we  obtain: 


l-yu2  +  v2  +  (uy  -  vx 


(.6  +  eh2  +  nk2)hkS>2 


)2]dxdy  s  y  Jj[u2+  y2  +  ^uy"Vx^'|dXdy 


s  :  e 


1(  It  f  D,  where  D  depends  on  the  region  only. 


JJl-yu2  +  v2  +  (uy  -  vx)2]dxdy  f  L2hk^2 


:■ „  , •  t-  <  f&We W. - 
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L  is  independent  of  h  and  k,  unless  they  are  of  the  same  order  for  any 

patch  of  rectangular  elements,  and  condition  (15)  is  proved,  since: 


sup 

(u,v)€UnxVh 


II 


’Kuy-vx)dxdy 


||(U.V)|| 


for  every 
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Appendix  B:  Oscillations  of  u  on  the  parabolic  line  using  numerical  projection 


1)  Pos s ible  reason  for  the  numerical  oscillations 


■a 


The  explicit  schemes  in  the  hyperbolic  region  with  Cauchy  boundary  conditidhs, 
using  functions  which  are  linear  in  x,y  ,  are: 


(la) 


(i,j) 


X  .  \tl  iJillLtl 

i.l 


T 


CiJ-l) 


IH 


Uc) 


Vi.j  =  -  0-5 


/i, j-l+vi+l,j-l]  +  jJXi+l,j-l  “  Xi,j-l) 


u.  .  =  0.5 


[3] 2/3 

ui.M*  ui*u-i)  *  “w  (j2/I-o-n,/s)(vi+1.j-rTt.j.i] 


u  and  u.  .  .  .  have  the  same  coefficient  and  the  same  sign  in  all  the  schemes. 

1  *  J  ‘  1  1 + 1  *  J  ~  1 

Therefore  u.  n-d  and  u.  .  n+d  lead  to  the  same  solution  in  the  hyperbolic  region. 

1+1, u 

The  same  consideration  cannot  be  done  for  v  and  \  because  they  have  different  signs 
in  some  of  the  schemes. 

Thp  numerical  projection  is  obtained  by  substituing  the  boundary  conditions 
along  Prefigure  2): 


“-N.j  ^  *  V-N > j  ' 


M 

!l 

i 

id 


_ ^  A.  ■ 


where  u  „  .  is  obtained  from  (lc)  and  v  „  .  from  (lb) . 

-N»j  -N*i 

3  2/3 

For  the  first  step,  j  ■  1,  and  y  =»  (^h)  ,  we  project  the  condition 

3  1/3 

u  N  1(3*0  +  v  N  1  *  f(y)  on  the  segment  {-Nh  $  x  s  (-N+l)h,y  ■  0  } 

The  values  of  u  N  Q  and  X_N  Q  are  known  (v_N  Q  ■  f(0)  and  u_N  is  computed 
by  L'  Hopital's  rule  and  X  is  taken  as  an  arbitrary  value  at  this  point),  therefore 
the  first  projection  gives  a  connection  between  u_N+1  Q,  V_N+1  0  and  X-n+1  0* 


(2) 


The  coefficient  of  u  .  n  decreases  to  zero  when  h  -+  0,  contrary  to  the 

-N  +  1 ,U 

coefficients  of  v  „  ,  _  and  X  .  ~  For  a  final  value  of  h  equation  (2) 

-N+1,U  -N  +  1,U* 

enables  a  large  error  o  u  n+1  0  an<*  error  may  increase  if  finer  meshes  are 
used. 


Assume  that 


u.  n  and  u.  .  _ 
1,0  1+1,0 


u  x.  ,  .  =  u^11?1!  +  d.  As  we  have  seen  before  we  can  add  ±d  to 
-N+1,0  -N+1,U 

accordingly,  without  changing  the  solution  in  the  hyperbolic  region. 


Therefore  oscillations  occure  in  the  solution  of  u  on  y  *  0. 
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2*  Explanation  why  good  results  are  achieved  with  A  piecewise  constant 
but  bad  convergence  using  A  bilinear  in  the  elliptic  rc2ion 
The  boundary  conditions  for  the  points  on  the  parabolic  line  are: 

i)  the  numerical  projection 

ii)  3J  (u,v,A)  *  0 
3vi,0 

iii)  3J  (u,v, A)  »  0 


where  J  is  the  functional  in  the  elliptic  region; 


C-N.O)  C-N+j.O) 


(-N,-l)l 


C-N+1,-1] 


.  (i-i.,Q)Cfc»oKiti.Q)  . 


:i-l,-l) (i,-l| (i+1,-1) 


(N,0) 


(-N.-N) 


(N,-N) 


The  values  of  u^  Q  appear  in  the  numerical  projection  and  in  the  schemes 

3J  =0  and  3d  =0  (and  3J.  if  ft  is  bilinear).  JW  =  0 
3ui,-l  3Xi,0  3Xi,-l  3ui,-l 

gives  the  following  coefficients: 


-h  <4'Y.O  *  ui+l ,0  *  * 


Therefore,  if  there  are  oscillations,  they  are  multiplied  by  frh  and 
do  not  influence  the  results  for  interior  points  for  small  h.  On 

the  other  hand,  ail  »  0  ^ 

a  Aj  n 


« 


(a)  X  bilinear 


4  (\,e  ‘  \-i>  *  tui .  1,0 '  “i-i.-i)  *  (ui-i,o  ‘ 


2  tvi  *  1,0  ‘  Ti-1,0)  *  tvi  ♦  1,-1  '  Vl,-1> 


(b)  X  piecewise  constant 


U  ij9  "  ui,-l  +  ui  +1,0  ui  +  1,-1  3  vi  +  1  vi,0  +  Vi  +  1,-1  ’  Vi,-1 


In  case  (a)  u.  ,  ,  u.  *  and  u.  .  n  have  different  coefficients. 
Therefore,  oscillations  on  the  parabolic  line  will  disturb  the  results 
in  the  elliptic  region. 

In  case  (b)  u  .  and  u.  .  have  the  same  coefficient;  hence,  the  oscil- 

X  p  ('  1  ♦ 

lations  along  the  parabolic  line  do  not  affect  the  other  variables  inside  the 
region. 

3 .  Numerical  results 

If  X  is  piecewise  constant,  the  oscillations  on  y  =  0  can  be 
decreased  by  fixing  the  arbitrary  value  of  X  at  the  point  (N,0) 
instead  of  (-N,0).  Then  the  equation  for  the  point  (-N,0)  is: 


8A-N,0 


*  0  ^  U-N,0  ‘  U-N,-l  +  U-N  +  1,0  U-N  +  1,-1 


=  v  —  v  +  v  —  v 

-N  +  1,0  -N,0  -N  +  1,-1  -N,  -1. 


We  know  the  values  of  u  _,  v  „,  v  -  from  the  boundary  conditions, 

■N 1 0  ” N | U  *N  |  “ 1 

so  we  have  the  following  connection: 


U-N  +  1,0  "  ll-N  +  1,-1  "u  -N,  -1 


-N  +  1,0 


-N  +  1,-1 


In  this  scheme  the  coefficient  of  u  _  does  not  decrease  to  zero 

-N  +  1,0 

when  h  -+  0  in  contrary  to  (2).  Therefore,  the  error  d  will  be  smaller. 


-  11  '!  'P'lP 
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Notes : 

1.  In  the  elliptic  region  with  Dirichlet  conditions  and  X  piecewise 

constant,  we  do  not  use  the  value  of  XN  Q  because  we  take  in  every 

rectangle  the  value  of  X  at  the  left  and  upper  corner.  But  if 

the  elliptic  problem  is  solved  with  the  numerical  projection,  then 

there  is  a  scheme,  belonging  to  the  hyperbolic  part,  which  contains 

n»  vm  n  ai*d  XKt  n.  In  this  case  we  fix  the  value  of  X„  If  X  „  A 
NjO  NjO  N,0  N | 0  “NjO 

is  given  as  an  arbitrary  value,  we  must  add  an  equation  for  X^  Q. 

We  take,  for  example,  extrapolation  of  X^  ^  q  anc*  ^  2  0‘  as 

we  have  explained  before,  larger  oscillations  appear  in  this  case, 

2.  If  X  is  bilinear  in  the  elliptic  region,  4  values  of  X  are  given 
around  the  corner  (-N,  -N)  in  both  cases;  the  Dirichlet  problem 

and  the  problem  with  the  numerical  projection. 

In  Table  VI  the  results  are  illustrated  for  the  first  example. 

In  this  case  u  ^  l*»0)  -  One  can  see  that  the  size  of  the 
oscillations  do  not  depend  on  h  but  they  increase  when  x  approaches  1. 
The  solution  does  not  converge  locally  on  y  ■  0  but  converges  in 


L 2  -  norm  and  in  the  weighted  norm  ( 

4.  Correction  on  the  parabolic  line 
We  compute 


2  1/2 

yu  dxdy)  '  (Table  V), 


di 


u.  .  -  2u.  +  u.  ,  . 

1-1,0  1,0  1  +  1,0 


Vi  -N<i<N 


where  Q  is  the  solution  on  y  *  0  of  the  elliptic  problem  with  the 
numerical  projection. 


-  Jk  \ 


I 

'I 
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We  define  new  boundary  conditions  on  the  parabolic  line: 

U-N,0  “  U-N,0  d-N  +  1 

"i.O  =  Ui*0  “  di 

Xo  ■  Xo  "  dN-l 
In  the  first  example  we  have: 


0.68 

< 

Kl 

< 

0.98 

h  - 

1/4 

0.49 

< 

Mil 

< 

1.03 

h  * 

1/8 

0.33 

< 

Mil 

< 

1.01 

h  - 

1/16 

and  the  smoothed  values  |u^  Q|  <_  0.06. 

We  solve  the  elliptic  region  again  but  with  u*  Q  as  Dirichlet 
boundary  conditions  on  y  ■  0.  XN  Q  does  not  appear  in  the  Dirichlet 
problem  if  X  is  piecewise  constant;  therefore,  1  1  ^ro,n 
previous  solution  is  taken  as  an  arbitrary  value  for  X.  The  results 
are  given  in  Table  V  , 


-  76  - 


Legend 


EU  - 

^  ^Unum  "  uanaP 

ixW1^ 

EUy  - 

(Iy(unum  “  uanal"* 

W)1'2 

EV  « 

^vnum  "  vanal^ 

LW)U2 

EG  - 

<E'X„u»  -  w2 

h  ■  1/4 

h  -  1/8 

h  «  1/16 

Rate  of 
Convergence 

EU 

0.0107 

0.0025 

0.0006 

2.15 

EV 

0.0453 

0.0104 

0.0025 

2.09 

EG 

0.2270 

0.0549 

0.0135 

2.05 

Table  la: 

Results  of  the  errors  and  rates  of  convergence  for  the 
Dirichlet  problem  using  X  bilinear. 

h  =  1/4 

h  =  1/8 

h  =  1/16 

Rate  of 
Convergence 

EU 

0.0629 

0.0257 

0.0100 

1.33 

EV 

0.0962 

0.0362 

0.0133 

1.43 

EG 

0.334 

0.184 

0.095 

0.91 

Table  lb:  Results  of  the  errors  and  rates  of  convergence  for  the 
Dirichlet  problem  using  piecewise  constant. 


h  *  1/4 
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h  =  1/16 


Rate  of 
Convergence 


h  »  1/8 


EU: 


Goursat 

0.116 

0.045 

0.017 

1.38 

Cauchy 

0.00289 

0.00092 

0.00033 

1.63 

EUy: 

Goursat 

0.091 

0.034 

0.012 

1.53 

Cauchy 

0.00269 

0.00086 

0,00028 

1.63 

EV: 

Goursat 

0.098 

0.035 

0.012 

1.52 

Cauchy 

0.0328 

0.0093 

0.0028 

1.77 

EG: 

Goursat 

0.006 

0.0022 

0.0006 

1.69 

Cauchy 

0.026 

0.0093 

0.0033 

1.49 

Table  II:  Results  of  the  errors  and  rates  cf  convergence  in  the 
hyperbolic  domain  with  exact  boundary  conditions  for 
ei  .her  the  Goursat  or  Cauchy  problem. 
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Rate  of 


H 

n 

JS 

h  =■  1/8 

h  «  1/16 

Convergence 

EU: 

y  <_  0 

0.265 

0.079 

0.020 

1.86 

y  >  0,  Goursat 

0.262 

0.083 

0.027 

1.64 

y  >  0,  Cauchy 

0.160 

0.111 

0.082 

0.48 

EUy: 

y  <_  0 

0.178 

0.059 

0.016 

1.73 

y  >  0,  Goursat 

0.212 

0.064 

0.020 

1.70 

y  >  0,  Cauchy 

0.138 

0.086 

0.059 

0.61 

EV: 

' 

y  <  0 

0.176 

0.057 

0.015 

1.88 

y  >  0,  Goursat 

0.269 

0.076 

0.022 

1.80 

y  >  0,  Cauchy 

0.169 

0.095 

0.062 

0.71 

EG: 

\  : 

y  <.  0 

0.634 

0.233 

0.067 

1.62 

y  >  0,  Goursat 

0.534 

0.202 

0.058 

1.60 

y  >  9,  Cauchy 

0.127 

0.048 

0.012 

1.70 

Table  Ilia:  Results  of  the  errors  and  rates  of  convergence  in  the 
mixed  domain  using  Bitsadze  relation  (11) . 

y  <_  0  X  bilinear 
y  >  0  Goursat  or  Cauchy 
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Rate  of 


h  ■  1/4 

h  =  1/8 

h  *  1/16 

Convergence 

EU:. 

y  <  0 

0.183 

0.066 

0.018 

1.67 

y  >  0,  Goursat 

0.339 

0.160 

0.054 

1.33 

y  >  0,  Cauchy 

0.176 

0.071 

0.024 

1.44 

EUy :  • 

y  <  0 

0.071 

0.027 

0.009 

1.49 

y  >0,  Goursat 

0.285 

0.129 

0.042 

5  .  34 

y  >  0,  Cauchy 

0.161 

0.064 

0.022 

1.44 

EV: 

y  <  0 

0.119 

0.042 

0.014 

1.55 

y  >  0,  Goursat 

0.335 

0.140 

0.044 

1.47 

y  >  0,  Cauchy 

0.239 

0.028 

1.55 

EG: 


y  <  0 

0,452 

1.07 

y  >  0,  Goursat 

0.361 

0.033 

1.23 

y  >  0,  Cauchy 

0.059 ■ 

0.008 

1.44 

Table  Illb:  Results  of  the  errors  and  rates  of  convergence  in  the 
mixed  domain  using  Bitsadze  relation  111). 

y  <  0:  \  piecewise  constant 

Goursat  or  Cauchy 


y  >  0 
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Rate  of 


h  =  1/4 

h  =  1/8 

h  «=  1/16 

Convergence 

EU: 

y  <_  0 

0.611 

0.426 

0.288 

0.54 

y  >  0 

0.224 

0.098 

0.041 

1.23 

EUy : 

y  <_  0 

0.081 

0.038 

0,018 

1.08 

V 

o 

0.172 

0.070 

0.028 

1.31 

EV: 

y  <  0 

0.143 

0.060 

0.026 

1.23 

y  >  0 

0.172 

0.070 

0.028 

1.31 

EG: 

y  <  0 

0.374 

0.174 

0.083 

1.08 

O 

A 

0.123 

0.038 

0.012 

1.68 

Table  IV:  Results  of  the  errors  and  rates  of  convergence  in  the 
mixed  domain  using  numerical  projection. 

y  <_  0:  A  piecewise  constant 

y  >  0:  Cauchy 
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Rate  of 


h  =  1/4 

h  =  1/8 

h  ■  1/16 

Convergence 

EU: 

y  i. 

0 

0.133 

0.071 

0.037 

0.92 

y  > 

0 

0.053 

0.029 

0.015 

0.91 

EUy : 

y  1 

0 

0.076 

0.C36 

0.016 

1.13 

y  > 

0 

0.050 

0.026 

0.012 

1.03 

EV: 

y  1 

0 

0.124 

0.052 

0.022 

1.25 

y  > 

0 

0.192 

0.068 

0.024 

1.50 

EG: 

y  1 

0 

0.262 

0.123 

0.061 

1.05 

y  > 

0 

0.035 

0.009 

0.002 

2.00 

Table  V:  Results  of  the  errors  and  rate  of  convergence  in  the 

mixed  domain,  usi.ng  numerical  projection  and  smoothing 
of  the  solution  on  the  parabolic  line. 

y  0:  A  piecewise  constant 

y  >  0  Cauchy 


■M 
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UNIQUENESS  OF  THE  SOLUTION  OF  THE  LAPLACE  EQUATION 
WITH  SPECIAL  NON-LINEAR  BOUNDARY  CONDITIONS 


Sara  Yanlv 


Abstract 

Non-standard,  non-linear  boundary  conditions  for  the  Laplace 
equation  are  considered,  and  the  question  of  uniqueness  examined. 
The  nonlinear  conditions  prescribe  the  mgnitude  of  the  gradient 
on  part  of  the  boundary,  and  with  either  Dirichlet  or  Newmann  con¬ 
ditions  prescribed  on  the  remaining  part.  The  solution  of  the 
problem  is  unique  up  to  a  sign. 
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^mA»ry  value  problem 


Consider  Laplace  equation: 


AU  =  O 


Cla) 


in  a  two  dimensional  domain  n  ,  bounded  by 


an  =  a^  u  an2 

with  the  non-standard  boundary  conditions: 


U|  =  0 
a»i 


|7U| 


1 


Ob) 


The  solution  for  problem  Cl)  is  not  unique;  the  question  is  how  many 
solutions  there  are  and  what  is  the  connection  among  them. 

Theorem  1: 

If  U^Cx,y)  is  a  solution  of  problem  Cla),  Clb),  then  U^Cx,y)  and 
-U^(x,y)  are  the  only  solutions  of  the  problem. 
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Proof: 


Using  complex  functions  of  a  complex  variable  z  =  x  +  iy  we  define  the 


function 


fCz)  =  4>(x,y)  +  i4)(x,y), 


ffe)  is  an  analytic  function  satisfying  the  following  boundary  conditions: 


Re  f  |  *  o 


Re  df'j  =  Re  f«  dz|  =  0 
30i  dfl! 


(which  is  equivalent  to  dU  =  l^dx  +  U  dy|  =  0) 

90x 

anH  I  -PI  I  -  1^1  -  l 


Let  us  test  the  problem  for  the  analytic  function  f’(z)  *  uCx,y) 
satisfying  conditons  (2*  )  and  (3)  , 


-  iv(x,y) 


Assuming  fj|  and  f£  are  analytic  functions  satisfying  (2')  and  (3),  then 


f’l  =  f£eiaCz) 


Im  a(z) J  =  0 

sn2 


;■&**»* 


■a Aiikio i  i. .iaiminMili.1  ,ii i,  n  .  ■■■-  id  i--niL-.it  Jiiaiim«yttaUfeiAn-iiai 


a(z)  =  a(x,y)  +  Ib(x,y) 


4tl.ni  .in,  u 
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Another  problem  is: 

AU  =  0  (5)  in  fl  which  is  bounded  by  3^i  U  882  with  the 
boundary  conditions: 

(5a) 

(5b) 

(x0,y0)  €  n 

Again  we  lave  the  same  theorem: 


Theorem  2: 

If  U  =  U  ^  is  a  solution  of  problem  ( 5 ) ,  then  U  =  are 

the  only  solutions  of  the  problem. 


and  along  38  ^  the  following  holds; 


so. 


since 


and 


Im  f„  dz|  =  Im  f,'  e^^dzj  =  ie  ^  sin  a(udx  +  vdy) 
1  30i  X  8Qi 

sin  a(x,y)  =0  ■*  a(x,y)  =  nir  . 

da(x,y)  =  axdx  +  aydy  =  bydx  -  lyiy  =  |£|  =0 

38! 

Kx,y)|  =  0 
..  3122 


we  get  b(x,y}  s  0 

and  a(x,y)  =  nir  in  S, 


hence 


UNIQUENESS  OF  THE  SOLUTION  OF  ELLIPTIC  PROBLEMS 
WITH  MIXED  BOUNDARY  CONDITIONS 


Sara  Yaniv 


Abstract 

In  the  problem  of  designing  shock-free  supercritical  air¬ 
foils  for  the  transonic  small  disturbance  theory*  certain 
mixed  boundary  conditions  arise. 

In  the  following,  the  question  of  uniqueness  of  a 
linear  elliptic  problem  with  the  appropriate  boundary 
conditions  is  examined. 


tfWWVi* 


1 ,  Introduction 

The  most  important  question  in  the  theory  of  elliptic  equations 
is  the  study  of  boundary  value  conditions  for  which  the  problem  is 
well  posed. 

The  most  investigated  problems  are  the  first (Dir ichlet) , second 
(Neumann  and  third(mixed)  boundary  value  problems.  All  these  problems 
are  special  cases  of  the  boundary  conditions: 


Y<rr  t6u*4> 


along  the  boundary 


see  [2]  . 

In  what  follows  we  prove  that  if  the  oblique  derivative  condition 

8-  Shi* 

along  the  boundary  is  given, the  Helmholtz  equation  is  unique  if 
a<0  .which  means  that  1  is  oriented  towards  the  exterior  of  the 
region  [2] . 


-L.  * 
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2 .  Boundary  Value  Problem 

Consider  the  self-adjoint  equation 

(1)  Au  *  Pu  in  ft 

for  P  =  P(x1,  x 2 1 • • • ,xn)  >  0. 

The  boundary  condition  is: 

In  +  a  I?  a  £  along  * 

where  is  the  outward  normal  derivative  and  g-g-  is  the 

tangential  derivative  along  the  boundary.  This  condition 
is  similar  to  the  known  mixed  boundary  condition  for 
elliptic  problem  [1]: 

(3)  —  +au  *  f  along  3ft. 

In  the  following  theorem  we  prove  that  problem  (1),  (2) 
has  a  unique  solution  under  the  condition  a  <  0,  (the 
same  holds  for  problem  (1),  (3)). 

Theorem: 

If  a (x^ ,X2 . . . ,xn)  <  0  then  problem  (1),  (2)  has  a  unique 


solution. 


Proof 


97  - 


Using  Green’s  identity  we  have: 


(4)  / 

fl 


/ 

3ft 


3u  3v 
3x^  5x ^ 


.  3u  3v  K 

+  Sx7  S3 r  +  uAv 
n  n 


dx,  . .  ,dx„ 
a  n 


r  n 


,3-1  „3v 


dx i**i dx  *  *  * dx  *  * • « dx _ 


u5xT  j - n 

«.  4  J 


Taking  v  a  solution  of  equation  (1)  and  v  »  u, 
then  (4)  reduces  to: 


(5)  l 


'dr-,2  »  (9U 

'“3X,J  ^ 


2 

Sxj  >  * 


c|S-i  *  p»2 


dx 


Let  us  introduce  the  notation: 


«II  -  / 

Si 


<iir>  *  f«2' 

l  n 


dx 


since  P  >  0  this  relation  can  be  identified  as  a  norm. 

Consider  two  solutions  and  u2  of  the  problem  (1),  (2)  in 
ft.  Suppose  u1 (x^ , . . ,xn)  and  u2(x1,...x  )  are  smooth  functions 
so  that  formula  (5)  can  be  applied  to  the  difference  u  *  u^  *  u2* 


T 
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Since  u  satisfies  (1)  in  8  we  have 
(6)  ||u||2.-9/u|HdS, 

and  from  equation  (2) 

.  -a  |2» 

3n  3s 

then 

null2 .  iou||  dS  -  \  fo  da. 

an  *  3n 

Using  the  mean-value  theorem  for/^  smooth  function  a (x^ , . . .  ,xn) , 
we  have : 

u  n 2  -  «  -  o. 

f  i2 

since  ds  =*  m  >  0 

J  3s 
3  Q 

then 

!  lu  1 1  2  -  ^  am  *  0 . 

Assuming  a  <  0  we  obtain: 


Mi  »  0 


and 


Thus,  uniqueness  is  assumed  for  problem  (1),  (2)  provided 
a  <  0  and  P  >  0. 

If  pC*i . v  *  0 

then 

U1  -  u2  -  c. 

which  means  that  under  the  boundary  condition  (2),  the 
solution  of  the  Laplace  equation 

Au  «  o 

can  be  determined  up  to  an  additive  content. 


We  solve  Laplace  equation  in  polar  coordinates: 

(6)  r2  $rr  ♦  r4»r  ♦  4>Qe  ■  f(r,6)  in  t) 

where 

ft  •  {(r,e)/rQ<  r  <  R,  0  «  6  < 

with  the  following  boundary  conditions: 

(7a)  ||-  (r,e  -  0)  -  gx(r) 

(7b)  $(r,$  -  J  )  -  g2(r) 

(7C)  ♦  (R.O  -  hx(9) 

(7d)  co*@-||5lSJ.  ■  hz(9) 

°  r  *  ro 

(figure  1) . 

Condition  (7d)  is  of  the  form  (2a)  for  which; 

«  •  -  0 
0 

0  <  9  <  7 


for 
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The  numerical  calculations  are  performed  for  the  following 
examples : 

I.  The  functions  f(r,6),  g^Cr),  g2(r),  hj  (6)  and  h2(9) 
are  chosen  in  such  a  way  that  the  function 

(8)  4>(r,6)  ■  r(cos8  *  sin  6) 

is  the  solution  of  the  problem. 

The  numerical  solution  converges  to  the  analytic  solution  (8) 
Error  estimates  and  numerical  rate  of  convergence  are  given  in 
table  (I) . 

II.  The  above  functions  were  chosen  so  that  the  function 

<J>(r,6)  *  1/r  is  the  solution  of  the  problem.  Results  are 


given  in  table  (II)* 
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Remark  ! 

Numerical  experiments  show  that  the  condition  a  5  0  is 
also  necessary  for  uniqueness. 

We  tried  to  solve  problems  I  and  II  for  the  half  disc: 

/ 

*  {(**6)  /rQ  (  r  (  R,  0  (  0  (  it} 

for  which 

a  - 

ro 

changes  sign  in  the  interval  0  *  6  <  it. 

In  this  case  the  solution  is  not  unique  and  we  are  unable 


to  get  a  solution. 


t 


Let 


A6  « 


IT 


Ar 


where 


R  *  1.6 


r  *  0.08 
o 


E  ,  » 

anal 


[s(4> 


2  I1'2 

anal  '^num^  ^®^rJ 


IDEG  - 


order  of  approximation  of  the  normal  derivative  on 
the  boundaries  (first  or  -  second  order). 


i 


Table  I:  ^anal  “  x  +  y  a  r(cos6  +  sin0) 


N 

M 

IDEG 

EANAL 

19 

4 

1 

0.169 

19 

8 

1 

0.055 

38 

4 

1 

0.138 

38 

8 

1 

0.050 

38 

16 

1 

0.022 

19 

4 

2 

0.035 

19 

8 

2 

0.007 

38 

4 

2 

0.026 

38 

8 

2 

0.0064 

38 

16 

2 

0.0016 

rH  |  $-• 


Eanal 

17.3 

9.86 

3.66 

3.24 

1.39 


11.60 

5.49 

2.12 

1.59 

0.50 
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NUMERICAL  ASPECTS  FOR  A  NEW  DESIGN  PROCEDURE  OF 


SUPERCRITICAL  WINGS 


Sara  Yaniv,  Frieda  Loinger  and  Nima  Geffen 


Abstract 

The  free  sonic  line  method  for  two-dimensional  transonic  small 
disturbances  equation  in  the  subsonic  region  is  described  and  results 
of  numerical  tests  are  given. 
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1.  Introduction 

A  method  for  designing  shock-free  transonic  configura¬ 
tions  has  been  developed  by  Fung,  Sobieczky  and  Seebass  [1]. 

The  flow  about  a  given  wing  is  solved  for  modified  flow 
equations,  using  a  fictitious  gas  in  part  of  the  field  to 
guarantee  elliptic  conditions  everywhere,  except  on  a  para¬ 
bolic  line,  i.e.,  a  sonic  line.  The  velocity  components 
along  this  sonic  line  serve  as  initial  conditions  for  exten¬ 
sion.  into  the  supersonic  pocket  terminating  in  a  zero  stream¬ 
line  curve,  which  modifies  a  part  of  the  original  airfoil  to 
give  a  shock-free  supersonic  pocket. 

Skipping  the  fictitious  gas  idea  and  calculations  alto¬ 
gether,  the  proposal  here  is  to  prescribe  the  shape  of  the 
sonic  line  and  get  an  elliptic  problem  in  the  subsonic  region 
with  new  non-standard  boundary  conditions.  The  solution  in 
the  elliptic  region  and  along  the  parabolic  (sonic)  line  pro¬ 
vides  initial  values  for  the  hyperbolic  pocket.  The  completion 
of  the  shockless  flow  field  is  carried  out  by  continuation  of 
the  field  into  the  hyperbollo  pooket. 

As  a  first  step,  we  treat  the  transonic  small  disturbance 
equations,  with  boundary  conditions  given  about  a  olroular 
sonic  line  and  a  "circular”  wing. 
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2 .  Transonic  small  disturbance  theory  and  formulation  of  tht 

boundary  *a«.ga  problem. 

The  transonic  small  disturbance  equations  are  derived  by  an 
asymptotic  procedure  applied  to  the  exact  equations  of  gas  dynamics 
for  the  freestream  mach  number,  close  to  1. 

The  small  parameter  of  the  expansion  procedure  is  the  air¬ 
foil  thickness  ratio  6,  and  the  flow  is  represented  as  small  per- 
trubation  of  a  uniform  stream.  The  limit  process  is  associated 
with  the  approximation  6  o,  K,x,  y  fixed  where  y  »  y 

and  K  s  (1  -  M»2 )/ <$2^ 3  ,a  transonic  similarity  parameter  [2J. 

The  basic  transonic  equation  is: 

[K<frx  -  (Y  +  1>$x2/23x  ♦  *yy  3  °- 

or  alternatively 

(1)  CK  -  (y  ♦  1HX3*XX  +  ♦yy  5  °- 

Equation  (1)  is  hyperbolic  for  $  >  K/(y  ♦  1),  elliptic 

for  <|»x  <  K/(y  +  1)  and  parabolic  for  $x  =  K/(y  +  1). 

for  flow  around  a  wing,  conditions  of  tangent  flow  have  to 
be  prescribed.  For  the  first  (and  second)  approximation  the  con¬ 
ditions  in  the  plane  of  the  wing,  y  *  0,  are: 


-  no 


$p<x,0)  * 


F‘(x) 

0 


|x|  <  1 
|xj  >  1 


where  the  body  shape  is  given  by 

y  =  5  F(x) 

(the  geometry  of  the  problem  and  the  boundary  conditions  are  shown 
in  f igure  1 ) . 


In  the  design  problem  a  wing^  for  which  shock-free  flow 
occurs t is  to  be  constructed.  In  our  procedure  a  symmetric  sonic 
line  F  is  chosen,  resulting  in  the  following  boundary  value- 
problem  for  the  elliptic  flow  outside  the  sonic  line: 

(2)  ck  -  (y  +  D^xa  $xx  ♦  *  0 

(3)  ♦.  (x.O)  »  Ir’(x)  C  <  |x|  <  1 


(4) 


*  (x,y)|r  *  K/(y  ♦  1) 


(figure  2) 
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The  mixed  type  condition  (4)  can  be  written  'is: 

*x|r  5  In  *(x>  *  H  b(x)|r  *  xrhj 

where  |^j-  is  the  derivative  in  the  outward  normal  direction  and 

is  the  tangential  derivative.  The  boundary  conditon  (4)  is 
the  one  we  examined  in  [33 >  uniqueness  of  the  solution  does  not 
exist  for  every  T. 

For  a  numerical  solution  to  the  problem*  boundary  conditions 
for  the  far  field  are  needed  to  replace  the  free  stream  condition 
a  =  U  (1,0)  at  infinity. 

Regarding  equation  (1)  as  an  elliptic  equation  with  a  non-linear 
right-hand  side  we  have: 

L[$3  =  K*xx  +  =  [<y  ♦  D/23  (u2)x 


where  u  =  4>  . 

x 

Applying  Green's  formula  for  L  in  n,  assuming 
continuous  with  continuous  first  derivatives,  we  have: 

<*U*J  -  *  L[*3)d£dn  = 

Taking  <K£,n>x,y)  as  the  fundamental  solution  for  the  whole  plane 


*9$  ,9tfi 

*Tn  -  +7n 
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Li£  =  <5tx-£)  iCy-n) 

<P  =  sC£-*,n-y) 

and  using  a  similar  procedure  as  in  [2],  for  the  problem  C2)t  C3) 
and  (4)  but  omitting  the  possibility  of  shocks,  we  obtain: 

<Kx,y)  =  //u2(<j^  ♦  Ki*c)dCdn  ♦ 

^f.Sn  ('1'*  **)  *  *fln  +  |||]dl- 

dfi  a  r1  U  r2  U  r3  U  r4  U  r5.  (Figure  3). 


where 


Along 


«1  u  r2  U  r3  “  r4 


and  along 


rs  ** 


-/» (H  *  §*>  - 


2x  1 


x+Ky  2tt/E_£ 


where 


n  s  n(£)  describes  the  sonic  line, 


substituting  the  conditions: 


Ini 

'rjUr* 


F’(U, 


9  $ 
3n 


3  4>  3  <p 

5n>  *  5f 


|i,  -  Mr*D 


and  assuming,  Ts  is  a  symmetric  function  of  x. 


/  (^+i^*)d  n  =  0 

r,H 


we  get : 


I  -i  .uliLI  i.l iM  «UH <41  i  -Unl.^^L  |  t  I  f.  ,  r  _  (  ^ 


of  the  unknowns  of  the  problem. 
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3.  Circular  sortie  line 

The  transonic  equation  for  the  perturbation  potential  is 


(K  -  <Y*l)*x)*xx  ♦  =  0. 


Substituting 


x  =  *- 


we  obtain  the  equation: 


!  J 


(7) 


* —  +  r  Hi  ±  — 

XX  yyy  yx  'xx 


Equation  (6)  is  solved  for  tne  elliptic  part  of  the  problem 
by  assuming  a  given  circular  sonic  line: 


r5 :  x2  +  y2  =  rQ2  <  1  . 


In  polar  coordinates; 


x  =  r  cos  6 
y  =  r  sin  0 , 


equation  (6)  becomes: 


(8) 


0rr  +  r^r  +  ^7  ^00 


lii*  cos2o 

ml rr 


t  21221  +  *  SiJ'lL  *  sinfe,  + 

Yr0  r  y0@  __2 


rr  r 


U 


condition  (9c): 


0  <  6  i  it  )  does  not  satisfy  the  uniqueness  condition: 


tge  >  0  [3] , 

r 


An  additional  requirement  on  the  solution  is  needed  to 
prescribe  a  well  posed  problem.  Consider  an  odd  function 

x  ,  satisfying: 


4>  Cx  ,y )  s  -  (j,  (-x  ,y)  . 


It  is  obvious  that  such  a  function  is  a  solution. 


The  new  probelm  is 


find  the  solution  of  equation  C8) 


satisfying  the  following  boundary  conditions: 


(10a) 


4>0  (r,  6  =0) 


rr 


rrt  <  r  <  1 
o 

1  <  r  <  R_ 


(10b)  <J>  (r  ,6=j) 


(10c) 


<j>r  .cos0 


.  sine 
*eT~ 


0  <  e  i  tt/2 


D  cos9 
2ttK  r 

r  =  Fee 

0  <  0  <  TT/2. 


(lOd) 
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Numerical  procedure 

The  non-linear  problem  is  tested  for  some  boundary  conditions 
for  which  the  analytic  solution  is  known.  Therefore  we  take  a  non- 
homogeneous  equation  which  is  solved  by  the  following  iterations: 

(ll)  (l  -0,<s(n-1))*~(n)  *  *^(n>  =  fCS.y) 

Equation  (11)  is  transformed  into  polar  coordinates  with  the 
following  boundarv  conditions.' 

4>x(r0,e)  =  <j)r  cose  -  <fre  SiSi  *  gl(e),  o  <  e  < 

o 

*CRw,0)  =  g2C0)  ,  0  <  e  < 

♦0<r,e=O)  =  g3(r)  ,  r0  <  r  < 

4>(r,e=  ~)  =  g4(r)  ,  rQ  <  r  «  R. 

The  domain  ft’  =  {(r,e),  rc  *  r  *  R®  ,0*8*^-} 
is  devided  into  rectangles  of  size  A6  .  A r 

Roo"ro 

where  A8  =  and  Ar  =  — g — 

For  an  interior  point  ((i,j),  i  =  j  =2,...,!.') 


central  differences  are  used  in  the  equation. 
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On  0=0,  where  <$>  $  is  given  as  a  boundary  condition,  we  use  either 
a  first  or  a  second  order  approximation. 

On  r  =  r  ,  where  6  =  <J>  cos9  s^-n-  is  given,  a  central 

o 

approximation  for  <}>Q  is  taken  and  4>  is  either  a  first  or  second 

o  r 

orde-  approximation,  according  to  the  order  chosen  for  the  boundary 

6=0. 

The  schemes  described  beforv  lead  to  a  matrix  which  is  tii- 
diagonal  in  blocks  . 

If  first  order  approximations  are  used  on  r=  ro  then  the  size  of 

all  the  blocks  are  the  same :  (M+l)x(M+l)  and  contain  the  coefficients 

of  the  schemes  belonging  to  one  row  in  the  mesh. 

If  second  order  approx imaixons  are  used  on  r  =  rQ  then  each 

point  on  r  =  r  is  connected  to  the  rows  r  +Ar  and  r  +  2Ar. 

*  o  o  o 

Therefore  the  order  of  one  diagonal  block  will  be  2(M+1)  x  2CM+1) 
and  the  size  of  the  others  is(M+l)  x(M+l)  only. 

The  system  of  oquations  is  solved  by  a  direct  method  of  elimination 
by  blocks  which  needs  storage  of  twice  the  size  of 'the  greatest 
block  . 

The  numerical  procedure  is  stopped  when 
(t(*;n)  -  *(n-1>)2)1/2 


<  e  . 
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Numerical  results*. 


Let 


49  =  W 


Ar  = 


R  -  r 

oo  O 

N 


where  R  «  1 . 6 


r  m  0.08 

o 


Eanal  = 


^anal^num*  M  Ar 


1/2 


IDEG  -  order  of  approximation  on  the  boundaries 
(first  of  second  order). 


The  numerical  procedure  is  tried  for  the  functions  <J>  s  6  and 
$  =  r  for  which  the  schemes  are  exact. 

The  process  converges  for  every  initial  guess  and  the  error  decreases 
to  zero  for  every  size  of  mesh.  If  a  increases  the  procedure  is 
slower,  i.e.  more  iterations  are  neeaed  to  achieve  the  same  error. 

For  4>«x+y  the  errors  are  illustrated  in  table  I  and  for 
<J>*  r  cos(e/2)  in  table  II.  In  both  cases  the  results  do  not 
depend  on  the  initial  guess.  In  the  last  example  oscillations  are 
observed  on  r  =  rg  but  they  decrease  when  finer  meshes  are  used  and 
vanish  for  larger  values  of  r^.  In  order  to  achieve  an  elliptic 
problem  for  the  non-linear  equation  (11)  a  is  limited  by  1  in  the 
first  example  and  a*2/rQ  (“0.57)  in  the  second  example. 
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Table  I : 

N 

<J>  =  x+y 

H 

=  r(cos0  +  sinO) 

IDES 

Eanal 

a 

19 

4 

1 

0.170 

19 

3 

1 

0.056 

38 

4 

1 

0.138 

38 

8 

1 

0.050 

Q-.fll 

19 

4 

2 

0.035 

19 

8 

2 

0.007 

19 

16 

2 

0.0017 

38 

4 

2 

0.026 

38 

8 

2 

0.006 

19 

8 

2 

0.012 

19 

16 

2 

0.003 

0.7 

38 

8 

2 

0  .  Oil 

Table  II:  $ 


r1/2  cos  (6/2) 
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N 


M 


IDEG  E  anal  a 


19 

19 

38 

76 

38 


4 

8 

8 

8 

16 


1 

1 

1 

1 

1 


0.107 

0.045 

0.019 

0.014 

0.010 


19 

19 

38 

76 

38 


4 

8 

8 

8 

16 


2  0.027 

2  0.012 

2  0.0026 

2  0.0007 

2  0.0002 


19 

38 

76 


8 

8 

8 


2  0.0051 

2  0.00076 

2  0.00015 


o.s- 


38 


16 


0.00061 
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